NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  Celifornia 


A  MICROCOMPUTER-BASED  NETWORK  OPTIMIZATION  PACKAGE 

by 

Richard  Henry  Duff 
September  1981 

Advisor:  G.  G.  Brown 

Approved  for  public  release,  distribution  unlimited 


~  0 1  i  2 


SECURITY  CLASSIFICATION  0*  Tml  >MI  fllai  8 ma  Kmi—4) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
before  COMPLETING  FORM 

1  IIMIit  numocn  1.  oovt  accession  no. 

A  \)'A  'UJ1 

i  RCCiRiCnT's  catalog  numocn 

rvf 

4  TITLE  <m*  tukllllm)  * 

A  Microcomputer-Based  Network  Optimization 

Package 

*•  TVPl  OP  MPORT  M  »en<oo  covEnco 

Master's  Thesis 

September  1981 

i-  FERFORMINO  ONO.  REAORT  NUMOCN 

7.  AUTmONi'O 

Richard  Henry  Duff 

S  Rf  RFORMINQ  OOOANlZ ATION  NAME  ANO  AOOOESS 

Naval  Postgraduate  School 

Monterey,  California  93940 

10.  NNOONAM  CLEMENT  NROjECT  TAM 
AREA  0  WORK  UNIT  NUMBERS 

II.  CONTROLUNO  OFFICE  NAME  ANO  AOOOESS 

Naval  Postgraduate  School 

Monterey,  California  93940 

11-  •f*OHT  OATS 

September  1981 

If  NUMBER  OF  4*0(1 

130 

II.  MONITORING  AOENCV  name  A  AOOOESSflf  mtfrml  Mat  Canmlll nf  OfllcaJ 

IS-  SECURITY  class.  1*1  mi,  rSOanj 

Unclassified 

ISA.  OCCL  ASSl  FlC  ATIQn/ OORNGRAOinG 
SCnCOuLC 

it.  otSTmwrioN  ir atement  tSi  out  ffawanj 


Approved  for  public  release,  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (ml  i»m  mSmttmml  «im4  In  Blmak  **.  If  HHmnmt  frmm 


I*  SUBRLEMCNT ARY  NOTES 


IS.  KCV  10*01  rCMMM  w  m«m  iM  If  an*  Or  *<••* 

Networks,  microcomputer,  optimization,  linear  programming,  nonlinear 
programming,  mixed  integer  programing,  minimum  cost  network  flow, 
mathematical  Drogramming  assignment  model,  transportation  model, 
transshiDment  model,  fixed  charge  network,  nonlinear  network. 

10.  ABSTRACT  fCmlHwa  an  rmmmtma  •!*•  If  *•«••••»  MM  IfnHIf  *T  **•«*  —  Saif 

.An  important  branch  of  mathematical  programming  is  concerned  with  optimization 
in  systems  described  by  networks.  This  paper  describes  an  integrated  suite 
of  advanced  tecnniques  for  dealing  with  minimum  cost  network  flow  formulations 
Written  in  Pascal  and  implemented  on  a  microcomputer  reoresentati ve  of  current 
small  computer  technology  (the  APPLE  II),  this  package  places  unprecedented 
modeling  versatility  and  solution  capability  on  the  analyst's  desktop.  Able 


DO  1473  COITION  OF  I  NOV  ••  IS  OBSOLETE 


S/N  0  101*014*  SSO I 


SECURITY  CLASSIFICATION  0*  this  BASE  flkM  Sm  BntmmmB) 


gpcrmn  j  amnni 


to  solve  small  to  medium  size  problems  (3000  arcs  or  less)  at  reasonable 
speeds,  programs  to  handle  capacitated  linear,  nonlinear  (convex  separable), 
mixed  integer  and  elastic  ranged  linear  models  in  addition  to  comprehensive 
control  and  data  management  routines  are  included.  Problem  size  and  solution 
speed  benchmarks  are  given  for  a  variety  of  models.  _ 


I  Accession  For 

1  NT  IS  GRAM 
!  DTI:'  TAB 
|  Unannounced 


Ava i : cl ; i 1 ty  Cedes 
A nnd/or 
Lst  Cpeoi.il 


t  Forrn  1473 
>1  01^2-014-6601 


MeuMV*  ev*MN^lCA^I*M  §f  T»l*  P  tali’ 


am* 


Approved  for  public  release,  distribution  unlimited 


A  Microcomputer-Based 
Network  Optimization  Package 


by 


Richard  Henry  Duff 
Major,  United  States  Marine  Corps 
B.S.,  Drexel  Institute  of  Technology,  1968 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  OPERATIONS  RESEARCH 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September,  1981 


Author: 
Approved  by: 


f4.  c.  i\ mA  J*L  HUf.  A4  j 

/  X 

^  yGgBALD'G.  BROWN 


Thesis  Advisor 
Second  Reader 


Ch ai rman.  Department  of  Operatsote"  Research 

4 j  fyj 

of  Information  and  Policy  Sciences 


DearT^orf 


ABSTRACT 


An  important  branch  of  mathematical  programing  is  concerned  with 
optimization  in  systems  described  by  networks.  This  paper  describes  an 
integrated  suite  of  advanced  techniques  for  dealing  with  minimum  cost 
network  flow  formulations.  Written  in  Pascal  and  implemented  on  a 
microcomputer  representative  of  current  small  computer  technology  (the 
APPLE  II),  this  package  places  unprecedented  modeling  versatility  and 
solution  capability  on  the  analyst's  desktop.  Able  to  solve  small  to 
medium  size  problems  (3000  arcs  or  less)  at  reasonable  speeds,  programs 
to  handle  capacitated  linear,  nonlinear  (convex  separable),  mixed 
integer  and  elastic  ranged  linear  models  in  addition  to  comprehensive 
control  and  data  management  routines  are  included.  Problem  size  and 
solution  speed  benchmarks  are  given  for  a  variety  of  models. 
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NOTATION 


Except  as  otherwise  indicated  in  the  text,  the  following  notational 
conventions  have  been  used  for  all  mathematical  expressions: 


Scalars 

Vectors 

Vector  Components 

Matrices 
Matrix  element 

Transpose 

Sequences 


Lower  case  Greek  and  Latin  letters  (e.g.,  a,  a). 
Lower  case  Latin  letters  with  a  bar  above  (a). 
Lower  case  Latin  letters  with  a  lower  case 
Latin  subscript  (a^). 

Upper  case  Latin  letters  (A). 

Lower  case  Latin  letters  with  lower  case 
Latin  subscripts  (a^). 

Superscript  T  (AT). 

Enclosed  in  braces  ({a^}). 
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I .  INTRODUCTION 


Mathematical  programming  can  be  defined  as  the  use  of  mathematical 
representations  (models)  to  plan  (program)  an  allocation  of  scarce 
resources  among  competing  activities  [Ref.  1].  An  important  branch 
of  mathematical  programming  deals  with  optimization  in  systems  described 
by  networks  or  collections  of  points  (nodes)  connected  by  links  (arcs). 
Such  models  arise  explicitly  in  a  variety  of  applications  and  include 
many  familiar  distribution  and  transportation  problems.  Bradley  [Ref.  2] 
suggests  that  network  models  have  been  so  widely  used  because:  (1)  they 
accurately  model  many  applications,  (2)  they  are  more  readily  accepted  by 
nonanalysts  than  other  models  (they  pictorially  resemble  the  physical 
process  being  modeled),  and  (3)  efficient  algorithms  are  available. 
Additionally,  many  models  can  easily  be  transformed  into  equivalent 
network  representations  by  direct  manipulation  (e.g.,  assignment  problems) 
or  exploitation  of  a  primal-dual  relationship  (e.g.,  critical  path 
problems).  The  object  of  all  these  formulations  is  usually  to  minimize 
the  cost  of  moving  a  single  commodity  through  the  network.  The  general 
form  of  this  problem,  the  network  programming  problem  (NPP),  is 

(NPP)  MINIMIZE  f(x)  (cost  function) 

s.t.  Ax  s  b  (node  flow-conservation  constraints) 

1  <  x  <  u  (arc  flow  bounds). 
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where  x  is  an  n-dimensional  vector 
(n  is  the  number  of  arcs), 
x.  represents  the  magnitude  of  the  flow  on  arc  i, 

A  is  an  m  by  n  (m  is  the  number  of  nodes) 
matrix  defined  as  follows: 

+1  if  arc  j  is  directed  away  from  node  i 
a, •  *  -1  if  arc  j  is  directed  toward  node  i 

J 

0  otherwise, 

b  is  an  m-dimensional  vector  of  flow  requirements  at 
each  node,  and 

f(x)  is  some  function  that  relates  cost  to  arc  flow. 

Although  NPP  can  be  solved  by  classical  general  purpose  constrained 
methods  (depending  on  the  exact  form  of  f(x)),  better  techniques  exist. 
Specialized  data  structures  employed  in  conjuction  with  modifications 
of  traditional  optimization  procedures  have  resulted  in  the  development 
of  extremely  efficient  algorithms  for  the  solution  of  NPP  [e.g..  Ref.  3]. 
This  work  has  been  motivated  in  part  by  the  intrinsic  usefulness  and  wide 
applicability  of  these  models  and  further  stimulated  by  the  increased 
availability  of  computational  devices.  Due  to  these  advances,  the 
operations  analyst  is  now  able  to  represent  larger  network  models,  and 
answer  questions  concerning  their  optimal  flow  easier  than  ever  before. 

As  interest  in  network  programming  continues  to  grow,  the  supply 
of  applications  software  will  undoubtedly  keep  pace.  This  has  been  the 
case  to  date;  however,  the  vast  majority  of  the  emerging  software  appears 
very  specific  in  nature  and  tailored  to  particular  classes  of  network 
models.  Certainly  there  are  many  superb  codes  available,  but  if,  for 


example,  the  capability  to  work  with  both  linear  and  nonlinear  cost 
functions  is  desired,  then  the  employment  of  two  separate  and  possibly 
quite  different  programs  is  required.  Even  a  relatively  mundane  task 
such  as  linearizing  a  nonlinear  model  to  obtain  a  feasible  starting 
solution  point  either  requires  incorporation  of  extra  code  or  an  exercise 
in  data  manipulation  and  program  linking.  Such  shortcomings  indicate 
that  unified  network  flow  programming  packages  able  to  cope  with  several 
different  types  of  models  would  be  quite  useful. 

Concurrent  with  these  advances  in  network  programming,  computer 
technology  has  experienced  breath-taking  progress.  Not  only  have  com¬ 
puters  improved  in  sheer  power,  but  they  have  become  increasingly 
compact  and  less  expensive  (relative  to  capability).  A  dramatic  example 
of  this  data-processing  revolution  is  the  evolution  of  desk-top  computers. 
These  small  computers,  the  so-called  microcomputers,  appeared  on  the 
scene  in  the  early  1970 's  and  today  can  provide  computer  power  comparable 
to  some  room-sized  machines  of  a  decade  ago  at  a  tiny  fraction  of  the 
cost.  More  importantly,  these  devices  are  becoming  so  common  that  their 
accessibility  to  technical  personnel  is  forecasted  to  be  virtually 
universal  in  the  near  future  [Ref.  4].  E.  M.  1.  Beale  [Ref.  5]  commented 
on  this  dissemination  of  computer  power  and  predicted  increased  use  of 
microcomputers  by  operations  analysts  to  solve  smaller  models  locally 
with  reliance  on  computer  centers  for  large  projects.  Tanenbaum  [Ref.  6] 
suggests  that  the  cost  of  small  computers  relative  to  communication 
expenses  now  makes  it  attractive  to  analyze  data  at  its  source  and  send 
only  summaries  back  to  large  computers  via  networking  arrangements. 


Even  though  the  hardware  is  available  and  its  potential  clearly 
recognized,  suprisingly  little  operations  analysis  software  of  reasonable 
complexity  and  generality  has  been  reported  for  microcomputers.  Aside 
from  a  few  decision  analysis  programs  [e.g..  Ref.  7,  8,  9],  effort  has 
been  mainly  concentrated  on  statistical  applications.  Morgenson 
[Ref.  10]  and  Isbell  [Ref.  11],  for  example,  have  developed  extensive 
data  analysis  packages  incorporating  sophisticated  techniques.  Addition¬ 
ally,  there  are  numerous  commercial  statistical  products  of  varying 
quality  and  capability  for  microcomputers . 

Surprisingly,  microcomputer-based  mathematical  programming  software 
is  still  virtually  non-existant,  despite  the  obvious  utility  of  such 
programs.  Undoubtedly,  there  are  many  ad  hoc,  rudimentary  implementa¬ 
tions  of  basic  methods;  however,  no  reference  to  any  work  of  consequence 
has  been  found  in  the  open  literature.  Feasibility  studies  [Ref.  12,  13] 
verify  the  desirability  of  optimization  packages  for  small  computers,  but 
only  explore  rather  elementary  network  algorithms.  These  studies  further 
suggest  that  network  flow  problems  are  likely  candidates  for  microcomputer 
solution  because  of  the  efficient  algorithms  at  hand.  Economic  justifi¬ 
cation  for  employing  the  microcomputer  for  optimization  purposes  has  not 
been  established  conclusively.  One  attempt  to  favorably  compare  such  use 
with  the  alternative  of  large,  general  purpose  computer  systems  [Ref.  13] 
has  been  severely  criticized  [Ref.  14]  and  the  issue  remains  undecided. 

The  research  reported  here  investigates  the  construction  of  unified 
network  flow  optimization  packages,  the  mathematical  programming  potential 
of  microcomputers,  and  (secondarily)  the  economic  feasibility  of  such 
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devices  with  respect  to  optimization  applications.  This  work  was  under¬ 
taken  with  the  following  (superficially  disparate)  goals: 

1.  Development  of  a  unified,  versatile  network-flow  optimization 
package  capable  of  handling  a  variety  of  models  and  utilizing 
"state-of-the-art"  algorithms. 

2.  Implementation  of  this  package  on  a  widely  available  micro¬ 
computer  to  explore  the  usefulness  of  smaller  computers  in  an 
optimization  context. 

Portions  of  this  work  were  presented  at  the  CORS/ORSA/TIMS  joint 
meeting  in  Toronto,  May,  1981.  At  that  international  meeting,  attendees 
expressed  surprise  that  such  an  integrated  network  optimization  package 
used  a  microcomputer  as  a  host. 
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II.  DESIGN  AND  IMPLEMENTATION  CONSIDERATIONS 


Microcomputer  implementation  of  large  programming  projects  requires 
careful  consideration  of  every  aspect  of  program  design  in  order  to  fully 
exploit  the  limited  resources  available.  This  section  briefly  discusses 
the  design  criteria  utilized  and  outlines  the  rationale  for  the  decisions 
that  shaped  the  package. 

A.  BACKGROUND 

Invariably  many  different  program  configurations  can  be  constructed 
to  accomplish  a  given  task,  however,  some  approaches  are  better  than 
others.  Program  effectiveness  is  a  well-studied  field  and  the  literature 
contains  many  characterizations  of  superior  software  designs.  Kreitzberg 
and  Shneiderman  [Ref.  15],  for  instance,  define  a  "good"  program  as  one 
that  is  correct  (provides  desired  results),  fast,  accurate,  hardware  inde¬ 
pendent,  efficient  with  storage,  and  easily  modified.  Thenson  [Ref.  16] 
and  Lientz  [Ref.  17]  further  suggest  that  operations  analysis  software 
should  place  heavy  emphasis  on  "user  impact"  or  ease  of  use.  These 
attributes,  commendable  as  they  may  be,  are  of  little  value  if  the 
program  cannot  solve  meaningful  problems  in  a  reasonable  amount  of  time. 
Meaningful,  in  this  context,  includes  not  only  the  intrinsic  usefulness 
and  generality  of  the  modeling  facilities  provided,  but  also  the  size  and 
complexity  of  the  representable  formulations.  Since  it  is  highly  unlikely 
that  a  program  can  be  optimal  in  every  respect,  compromises  are  inevitable. 
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Numerous  factors  interact  to  establish  program  limitations  with 
the  inherent  capabilities  of  the  host  computer  playing  the  dominant 
role.  Choice  of  algorithms,  data  structures,  programming  language, 
memory  and  peripheral  management,  and  programming  tactics  are  also 
significant  determinants  of  ultimate  software  efficiency.  For  the 
traditional  user  of  large  computer  systems,  the  improper  choice  of  one  or 
more  of  these  variables  usually  results  in  merely  a  slightly  degraded 
package  rather  than  outright  fail  — e.  "’he  smal 1 -computer  environment  is 
less  forgiving  and  does  not  allow  ns  luxury  without  severe  performance 
penalties. 

B.  HARDWARE 

1.  Microcomputer  Characteristics 

A  typical  microcomputer  system  consists  of  a  central  processing 
unit ,  memory,  peripheral  devices,  and  software.  Exclusive  of  peripherals, 
such  packages  physically  resemble  an  electric  typewriter.  The  specific 
form  of  the  processing  unit  varies,  but  most  can  directly  address  65,536 
eight-bit  words  of  random-access  memory.  This  "fast"  (typically  200-400 
nsec)  or  core  memory  is  usually  supplemented  by  slower  (but  more  plentiful) 
storage  (worst  case,  1-2  seconds  seek  time)  in  the  form  of  disk  drives.  Each 
disk  drive  provides  approximately  one  hundred  thousand  (or  more)  words  of 
memory  space  on  removable  magnetic  media.  The  number  of  disk  drives  allowed, 
their  individual  capacities,  and  access  times  depend  on  the  microcomputer. 
More  sophisticated  (and  expensive)  offline  memory  devices  are  also  available, 
capable  of  storing  millions  of  words.  Common  peripherals  include  communica¬ 
tion  equipment  to  access  other  computers,  video  displays,  and  printers. 


2.  Target  Microcomputer 


Clearly,  as  technology  is  changing  so  rapidly,  any  attempt  to 
identify  an  optimal  hardware  selection  would  be  futile.  Instead,  an 
acceptable  system  representative  of  available  products,  the  APPLE  II* 
microcomputer,  is  considered.  All  software  development  has  been  per¬ 
formed  on  an  APPLE  II  system  configured  as  shown  in  Figure  1. 


1.  APPLE  II  microcomputer  with  65,000  words 
of  memory. 

2.  Two  disk-drive  offline  storage  devices 
(150,000  words  each). 

3.  UCSD  Pascal  (language  card). 

4.  Comnuni cat ions  device  (modem). 

5.  Printer. 

6.  Eighty  character  by  twenty-four  line 
video  display. 


Fig.  1.  Development  Hardware 

Although  this  machine  is  not  the  most  powerful  of  its  class,  it 
is  a  reasonable  choice  for  many  operation  analysis  applications  due  to 
its  availability,  capability,  and  low  cost.  Introduced  in  1977,  the 
APPLE  II  is  a  widely-used  device,  still  in  production,  and  likely  to 
remain  so  for  the  next  few  years  [Ref.  18].  As  over  200,000  units  have 
been  manufactured  and  distributed  worldwide,  a  broad  base  of  technical 
support  exists  and  numerous  firms  offer  peripheral  products.  The  APPLE 
II  can  support  all  of  the  popular  microcomputer  programming  languages 
(BASIC,  Pascal,  FORTRAN,  and  Assembly  language).  Easily  expandable, 

*APPLE  II  is  a  trademark  of  Apple  Computer  Inc. 
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complete  systems  can  be  purchased  for  between  two  and  four  thousand 
dollars . 

3.  Machine  Independence 

All  computers,  regardless  of  size,  have  unique  features  and  one 
must  be  careful  if  machine  independence  of  software  is  desired.  Given  the 
variety  of  microcomputers  available  today  and  the  rapid  changes  expected, 
portability  between  existing  and  conceptual  machines  is  necessary  if  the 
considerable  investment  involved  in  the  development  of  sophisticated 
programs  is  to  be  protected.  Thus  the  importance  of  machine  independence 
should  not  be  underestimated.  For  this  reason,  features  peculiar  to  the 
APPLE  II  have  been  avoided  wherever  possible.  Thus  the  software  presented 
should  run,  with  minimal  modifications,  on  any  equivalent  hardware 
package. 

C.  MODEL  AND  SOLUTION  TECHNIQUE  SELECTION 

The  determination  of  specific  processes  to  model  and  subsequent 
selection  of  solution  techniques  compatible  with  computing  resources  are 
critical  issues  in  the  design  of  any  optimization  software. 

1.  Models 

Once  the  general  class  of  problems  has  been  established,  a  few 
carefully  chosen  models  often  can  adequately  represent  the  majority  of 
anticipated  situations.  Specialized  models  can  then  be  added  as  necessary. 
The  following  minimum  cost  network  flow  models  are  considered  essential 
for  any  comprehensive  package: 

o  Linear  cost  function  with  bounded  variables, 
o  Nonlinear  convex  separable  cost  function  with  bounded  variables. 
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As  useful  as  these  two  models  are,  they  cannot  depict  some 
desired  formulations.  Additional  flexibility  is  provided  by: 

o  Linear  cost  function  with  bounded  variables  and  elastic  ranged 
constraints  (see  Section  III,  Subsection  C). 
o  Linear  cost  function  with  elastic  range  constraints  and  bounded 
variables,  any  of  which  may  also  be  specified  as  "1-u"  variables 
(only  allowable  flow  is  at  one  of  the  bounds). 

These  four  models  comprise  a  basic  package  capable  of  repre¬ 
senting  a  wide  variety  of  single  commodity  network  flow  problems. 

Included  are  fundamental  examples  of  linear,  nonlinear,  and  mixed-integer 
network  optimization  models.  Building  from  this  network  paradigm,  these 
features  can  be  used  to  represent  myriad  mathematical  models.  Other 
models  could  easily  be  added,  however,  time  constraints  precluded  further 
extension  of  the  package. 

2.  Solution  Methods 

When  dealing  with  microcomputers,  algorithms  are  difficult 
to  select.  Simple  algorithms  require  little  storage,  but  are  typically 
inefficient  and  slow.  More  complicated  approaches  offer  speed  at  the 
expense  of  increased  memory  usage  and  expanded  data  structure  requirements 
Large  computers  are  able  to  use  the  efficient  algorithms  because  the 
additional  storage  requirements  are  infinitesimal  compared  to  the  total 
memory  available.  For  small  machines,  the  choice  is  not  so  clear.  A 
large  portion  of  memory  is  consumed  by  the  added  complexity  of  the 
advanced  algorithms,  resulting  in  significant  reduction  of  maximum  problem 
size.  If  one  chooses  the  simpler  methods,  larger  programs  are  accommo¬ 
dated  but  processing  time  can  soar  to  unacceptable  levels  for  all  problems 


Execution  speed  appears  to  be  far  more  important  than  maximum 
problem  size.  Providing  quick  answers  to  small  problems  is  the  most 
likely  optimization  role  for  microcomputers  at  the  present  time.  Very 
large  problems  will  continue  to  require  extensive  computational  facilities 
well  beyond  the  capabilities  of  current  microcomputers.  Given  the 
predicted  advances  in  small  computers,  memory-saving  techniques  should 
decrease  in  importance  as  devices  with  greater  storage  capacities  become 
available.  Using  inefficient  methods  instead  of  the  more  capable  but 
memory- intensive  algorithms  does  not  seem  to  be  justified. 

Numerical  representation  must  also  be  considered.  Small  computers 
are  usually  more  restricted  than  their  larger  counterparts  in  (1)  mantissa 
precision,  (2)  real  variable  exponent  range,  and  (3)  maximum  allowable 
integer  size.  If  provisions  exist  (or  can  be  created)  to  deal  with  these 
problems,  performance  degradation  may  result.  Thus  candidate  algorithms 
for  microcomputer  use  must  be  relatively  insensitive  to  such  limitations. 

The  programs  selected  for  these  models  are  considered  state-of- 
the-art  and  all  represent  current  research  efforts.  Utilizing  advanced 
algorithms  and  efficient  data  structures,  these  programs  have  been 
included  on  the  basis  of  their  suitability  for  microcomputer  adaptation 
and  mutual  compatability.  The  original  programs  were  coded  in  FORTRAN 
for  implementation  on  large  computers.  Each  program  has  been  translated 
into  UCSO  Pascal  and  modified  as  necessary  to  enhance  its  efficiency  in  a 
microcomputer  environment. 

GNET  [Ref.  3]  is  the  foundation  of  the  package.  A  highly  regarded 
and  widely-used  code,  it  is  one  of  the  fastest  methods  available  for 
solving  linear  minimum  cost  network  models  with  bounded  variables. 
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Additionally,  GNET  can  be  used  to  linearize  nonlinear  problems  and 
produce  initial  feasible  solutions.  The  documented  speed  of  the  method 
[Ref.  2]  and  its  efficient  data  structures  make  it  ideal  for  microcomputer 
implementation. 

Nonlinear  formulations  are  solved  by  NLPNET  [Ref.  19],  a  program 
by  Dembo  that  utilizes  data  structures  in  the  spirit  of  GNET  and  new 
adaptable  direction-finding  techniques.  The  user  of  NLPNET  can  control 
the  accuracy  of  the  direction-finding  process  and  thus  the  amount  of 
computation  effort  expended  for  a  solution.  This  feature  makes  such  an 
algorithm  perfect  for  microcomputers. 

Elastic  ranged  constraints  are  handled  by  ENET,  a  recently 
introduced  code  by  Brown  and  Graves  [Ref.  20].  An  extension  of  GNET,  it 
uses  the  same  efficient  methods  and  data  structures. 

Finally,  UNET,  another  new  code  by  Brown  and  Graves  [Ref.  20],  is 
used  to  solve  elastic  ranged  constraint  problems  with  bounded  variables, 
some  of  which  may  be  specified  as  "1-u"  variables  (flow  at  either  bound, 
only).  Also  employing  data  structures  similar  to  GNET,  this  algorithm 
requires  successive  calls  to  ENET  to  solve  enumeration  subproblems. 

Although  ENET  can  duplicate  the  capabilities  of  GNET,  the  latter 
is  explicitly  retained  to  solve  larger  problems  of  the  simpler  structure 
amenable  to  GNET.  This  is  not  an  algorithmic  disadvantage  of  ENET,  but  a 
consequence  of  the  package  design.  Since  UNET  repeatedly  calls  ENET, 
both  procedures  are  kept  in  memory  together  to  reduce  disk  access  time. 
Therefore,  a<  the  programs  are  constructed,  GNET  is  a  more  compact  code 
than  the  ENET/UNET  combination.  It  is  intended  that  future  revisions  of 
the  package  will  employ  only  ENET  and  UNET. 
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D.  PROGRAMMING  LANGUAGE 


Pascal,  specifically  the  University  of  California  at  San  Diego  (UCSD) 
version,  appears  to  currently  offer  the  best  programming  language  facil¬ 
ities  for  large-scale  microcomputer  programming  projects.  Possible 
alternatives  are  BASIC  and  FORTRAN,  languages  that  suffer  from  crippling 
shortcomings  as  presently  implemented  on  microcomputers. 

1.  BASIC 

Once  the  only  high-level  language  available  for  small  computers, 
BASIC  remains  very  popular.  While  BASIC  is  perhaps  suitable  for  some 
applications,  serious  mathematical  programming  efforts  are  hampered  by 
its  limited  speed,  primitive  programing  power,  and  lack  of  transport¬ 
ability.  When  installed  as  an  interpreted  language  (the  usual  case), 
each  source  statement  must  be  translated  into  machine  language  every  time 
it  is  encountered  in  the  logical  flow  of  the  program.  This  results  in 
long  execution  times,  especially  for  complex  programs  with  many  iterative, 
structures  (i.e.,  optimization  programs).  Additionally,  most  dialects  of 
BASIC  allow  only  global  (accessable  to  all  program  portions)  variables 
with  restrictive  naming  conventions.  Variable  name  conflicts  often  arise 
in  such  circumstances  greatly  reducing  the  portability  of  subprograms. 

The  use  of  subprograms  is  further  hindered  by  the  inability  to  pass 
parameters  between  program  fragments  without  resorting  to  global  variable 
reassignments.  Compiled  versions  (one-time  translations  of  source 
statements  to  machine  language)  of  BASIC  are  available  that  supposedly 
alleviate  these  problems,  however,  they  are  necessarily  machine-specific, 
further  aggravating  the  transportability  issue.  There  are  probably  as 
many  different  versions  of  BASIC  on  the  market  as  there  are  types  of 
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microcomputers.  Hardware  manufacturers  modify  the  language  to  fit  their 
particular  needs  resulting  in  potential  difficulties  when  programs  are 
moved  to  different  machines.  It  is  indeed  a  rare  BASIC  program  that  can 
be  transported  without  modification,  with  the  chances  of  such  success 
diminishing  as  program  length  and  complexity  increases. 

2.  FORTRAN 

Implementations  of  FORTRAN  for  microcomputers,  although  increasing 
in  number,  are  still  rather  uncommon.  A  compiled  language,  FORTRAN 
compares  favorably  to  other  languages  in  terms  of  speed  and  programming 
power.  The  biggest  drawback  of  FORTRAN  is  lack  of  standardization,  a 
hinderance  that  has  plagued  FORTRAN  users  on  large  computers  for  years. 
Compiled  code  is  completely  machine-specific  and  thus  not  portable. 
Transportability  of  the  source  statements  varies  depending  of  the  degree 
of  similarity  between  the  two  versions  of  FORTRAN  in  question.  This  lack 
of  standardization  combined  with  the  inability  to  exchange  compiled  code 
between  different  microcomputers  makes  FORTRAN  a  poor  choice  if  program 
portability  is  a  goal . 

3.  Pascal 

a.  Standard  Pascal 

Pascal  is  a  relatively  new  language  (compared  to  BASIC 
and  FORTRAN)  having  been  formally  defined  by  N.  Wirth  in  1971  [Ref.  21]. 
Named  after  the  famous  mathematician  Blaise  Pascal,  the  language  was 
originally  intended  as  a  vehicle  to  teach  computer  programming.  An 
excellent  instructional  tool,  Pascal  is  now  recognized  as  a  powerful 
general-purpose  language.  Similar  to  ALGOL,  Pascal  provides  a  rich 
set  of  program  and  data  structuring  tools.  Thus  a  great  portion  of 
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b.  UCSD  Pascal 

These  impressive  programming  features  alone  are  enough  to 

warrant  the  use  of  Pascal  for  difficult  microcomputer  programming  projects. 

When  one  examines  the  UCSD  version,  however,  the  evidence  seems  over- 
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whelming.  An  extension  of  standard  Pascal,  UCSD  Pascal  was  developed 
by  a  group  headed  by  K.  L.  Bowles  in  1975.  Explicitly  designed  to  run  in 
a  small  computer  environment,  it  retains  most  of  the  features  of  standard 
Pascal  while  adding  many  extensions  that  increase  the  power  of  the 
language.  For  our  applications,  we  are  interested  in  two  of  these 
additional  facilities  that  (1)  enhance  memory  manipulation  capabilities 
and  (2)  provide  for  extremely  transportable  programs. 

The  memory  management  techniques  offered  by  UCSD  Pascal  are 
truly  remarkable  for  a  microcomputer  software  system.  The  programmer  has 
complete  control  over  which  portions  of  code  are  to  be  in  memory  at  a 
given  point  in  the  program  execution.  For  small  programs,  the  entire 
code  is  generally  loaded  into  memory  at  the  beginning  of  execution  while 
larger  programs  can  be  segmented  to  provide  the  desired  partitions.  A 
small  portion  of  the  program,  again  totally  definable  by  the  programmer, 
must  be  in  residence  at  all  times  to  control  the  overlaying  process. 
Complete  programs  can  be  called  into  memory,  each  with  its  own  local 
memory  control,  executed,  and  other  programs  subsequently  activated.  In 
this  manner,  programs  can  be  chained  together  to  form  coherent  packages. 
Additionally,  groups  of  often-used  blocks  of  code  can  be  pre-compiled  and 
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UCSD  Pascal  is  a  trademark  of  the  Regents  of  the  University  of 
California. 
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placed  in  libraries.  When  activated,  the  specific  library  code  is  loaded 
into  memory  with  the  calling  program  and  executed.  Library  routines  can 
either  remain  in  memory  for  the  entire  program  execution  process  or  be 
discarded  after  each  use  (at  the  programmer's  option).  These  methods 
allow  programs  that  would  otherwise  be  impossible  to  implement  on  a 
microcomputer--the  package  presented  is  such  a  case  as  the  combined 
code  of  the  component  programs  far  exceeds  the  memory  capacity  of  most 
microcomputers. 

UCSO  Pascal  is  a  hybrid  compi led/ interpreted  language,  a 
concept  that  produces  unprecedented  machine  independence.  Source  state¬ 
ments  are  first  compiled  to  an  intermediate  pseudo  (p)  code  that  is  then 
executed  by  a  machine-dependent  interpreter.  This  greatly  enhances 
portability  as  the  interpreter,  which  normally  takes  about  six  man-months 
to  write  [Ref.  22],  is  the  only  part  of  the  system  that  must  be  changed 
to  take  advantage  of  a  new  hardware  configuration.  The  original  compiled 
source  code  can  be  transported  without  modification  to  any  machine  that 
has  a  p-code  interpreter  installed.  Figure  3  depicts  the  relationship 
between  the  components  of  the  p-code  system  (Lewis  [Ref.  23]).  Such 
interpreters  now  exist  for  all  the  commonly  used  microcomputer  central 
processing  units  [Ref.  24].  The  use  of  p-code  is  not  restricted  to  UCSD 
Pascal;  other  languages  (e.g.,  FORTRAN-77  and  PL/ I)  have  been  treated  in 
this  manner  to  achieve  similar  portability.  However,  UCSO  Pascal  is 
currently  the  only  p-code  implementation  that  is  widely  available. 

As  previously  indicated,  speed  is  an  important  programming 
language  attribute.  The  hybrid  nature  of  UCSD  Pascal  results  in 
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Fig.  3.  Relationship  Between  P-Code  and  Pascal 

execution  time  somewhat  slower  than  pure-compiled  languages  such  as  FORTRAN, 
but  much  faster  than  interpreted  languages.  Gagne  [Ref.  25]  reports  that 
Pascal  runs  three  times  as  fast  as  the  very  best  BASIC  and  ten  to  thirty 
times  faster  than  most.  Informal  timing  tests  comparing  Pascal  with  APPLE  II 
BASIC  confirms  these  performance  assertions.  Any  loss  of  speed  through  the 
interpretation  of  p-code  is  more  than  offset  by  the  other  desirable  features 
of  the  language,  especially  portability, 
c.  Limitations  of  Pascal 

Pascal  is  not  without  its  difficulties--some  are  minor 
irritants  while  others  are  more  important.  Foremost  is  the  precision  of 
real  variables  allowed  by  the  APPLE  II  version  of  UCSD  Pascal.  Without 
taking  any  special  programing  measures,  six  to  seven  significant  figures 
can  be  accommodated  in  a  mantissa  representation  equivalent  to  IBM  360/370 
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single  precision.  It  must  be  understood  that  this  is  an  implementation 
issue  on  the  APPLE  II  and  not  an  intrinsic  limitation  of  Pascal.  By 
suitable  programming,  the  precision  can  be  extended;  however,  it  was 
decided  not  to  do  so  because  of  the  speed  degradation  that  would  surely 
have  resulted.  The  only  package  module  that  utilizes  vast  amounts  of 
real  arithmetic  is  NLPNET,  the  nonlinear  network  code.  As  shall  be 
shown,  for  the  nonlinear  problems  that  have  been  attempted,  limited 
precision  does  not  appear  to  hinder  performance. 

Another  disturbing  omission  is  that  all  arrays  are  static 
and  created  with  a  size  and  type  dope  vector.  Linkage  conventions 
require  that  actual  parameters  in  subprogram  calls  and  formal  parameters 
in  the  routine  called  must  agree  in  type  and  dimension.  This  precludes 
problem-dependent  memory  management  as  an  automatic  feature.  Since  array 
bounds  must  be  predeclared,  their  alteration,  for  any  reason,  requires 
complete  recompilation. 

Also,  compilation  itself  is  a  rather  slow  process.  The 
nominal  compilation  rate  for  APPLE  II  UCSD  Pascal  is  two  hundred  source 
statements  per  minute.  This  becomes  tedious  for  the  long  programs  devel¬ 
oped  for  this  package,  some  of  which  are  on  the  order  of  three  thousand 
source  lines  (exclusive  of  comments). 

The  final  important  criticism  is  the  lack  of  a  predefined 
exponentiation  operator,  a  crucial  primitive  for  nonlinear  optimization 
work.  For  development  purposes,  an  exponentiation  operator  has  been 
written  in  UCSD  Pascal.  A  much  better  solution  would  be  a  fast  machine 
language  function,  or  utilization  of  special  arithmetic  processing 
hardware. 
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E.  MEMORY  MANAGEMENT 

A  major  challenge  to  the  microcomputer  programmer  attempting  to  code 
a  complicated  algorithm  is  memory  management.  The  interaction  between 
the  two  types  of  memory  (fast  "core"  and  slower  peripherial)  is  a  deli¬ 
cate  affair  and  an  extremely  important  influence  on  overall  program 
efficiency.  As  a  general  rule,  optimization  programs  will  always  expand 
to  fill  available  core  memory.  Memory  not  occupied  by  program  code  will 
contain  the  data  structure  representation.  In  order  to  extend  problem 
size  further,  either  the  data  or  the  program  code  (perhaps  both)  must  be 
partitioned  in  some  manner  and  the  resulting  pieces  moved  in  and  out  of 
fast  memory  as  needed.  Both  these  approaches  are  feasible  with  today's 
microcomputers,  but  the  intrinsic  code  management  features  of  UCSD  Pascal 
make  program  code  segmenting  much  easier  to  implement. 

Data  partitioning  (so-called  "in-core/out-of -core"  operations) 
requires  additional  control  logic  and  data  structures  to  work  effectively. 
As  the  purpose  of  this  research  is  demonstrative  in  nature,  data  parti¬ 
tioning  experimentation  has  been  deferred  for  future  package  enhancement. 
Only  program  code  segmentation  has  been  utilized  to  dynamically  manage 
memory  resources.  Such  a  process  is  termed  "swapping"  or  "overlaying." 

Due  to  the  iterative  nature  of  mathematical  programming  algorithms 
and  the  time  expense  involved  in  accessing  current  disk  drives,  one  must 
be  extremely  mindful  of  how  the  program  is  segmented.  It  is  desirable  to 
have  minimal  program  code  in  residence  at  all  times  during  the  critical 
solution  process.  If  the  divisions  are  too  fine,  then  execution  speed 
suffers  as  the  disk  is  constantly  being  accessed.  Keeping  too  much  of 
the  program  in  core  reduces  disk  access  at  the  expense  of  data  storage. 
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Therefore,  each  program  must  be  examined  carefully  and  partitioned 
to  yield  the  best  compromise  of  speed  and  storage  management. 


F.  THE  USER  INTERFACE 

All  computer  programs  require  a  certain  level  of  human  interaction 
to  produce  results.  Depending  on  the  particular  application,  this 
participation  may  range  from  minimal  operations  such  as  data  input  and 
output  interpretation  to  complex  control  decisions.  Whatever  the  degree 
of  interaction  necessary,  software  must  be  designed  with  the  user  in  mind 
if  maximum  benefit  is  to  be  derived  from  its  use.  Programs  that  properly 
take  the  human  element  into  consideration  and  attempt  to  promote  meaning¬ 
ful  man-machine  dialog  are  termed  "user-friendly."  As  microcomputers  are 
interactive  devices  usually  operated  by  a  single  person  with  the  entire 
computing  system  at  arm's  reach,  user  involvement  is  inescapable. 

Software  designed  for  the  small  computer  must  therefore  be  user-friendly. 

1.  Data  Input 

Thenson  [Ref.  16]  states  that  programs  should  be  designed  to 
minimize  the  probability  of  errors  in  the  user  input  process.  Likely 
causes  of  these  errors  are  incompatible  input  format,  invalid  characters 
in  the  input  field,  and  inadmissible  input  values. 

Input  formats  should  be  flexible  whenever  possible,  allowing 
user  control  to  effect  format  reconfiguration.  Otherwise,  as  a  default, 
the  easiest  method  for  the  user  should  prevail  and  internal  conversion 
can  then  be  undertaken  as  needed. 

The  frustrating  problem  of  invalid  characters  in  the  input  stream 
can  be  avoided  by  the  use  of  buffering  and  filtering  techniques  to 
isolate  undesired  values.  This  error  is  often  encountered  during 
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keyboard  input  where,  for  example,  an  alphabetic  character  is  entered 
when  a  numeric  digit  is  required.  Most  operating  systems  respond  to  this 
difficulty  by  terminating  program  execution,  an  intolerable  situation. 
Morgenson  [Ref.  10]  describes  a  procedure  in  detail  where  all  numeric 
values  are  input  as  string  variables.  A  lexicographic  scan  on  the  string 
is  then  performed  to  identify  offending  characters  and  construct  the 
number  for  program  use.  We  use  a  modification  of  his  technique. 

Inadmissible  values  can  be  isolated  by  screening  all  input 
prior  to  its  use.  Preferably,  this  is  performed  as  the  data  is  received 
by  the  program.  Simple  precautions  such  as  confirming  that  the  data 
conforms  to  required  numerical  ranges  (range  screening)  and  sign  restric¬ 
tions  can  prevent  unexplained  abnormal  program  terminations. 

Every  attempt  should  be  made  to  not  only  detect  errors,  but  also 
to  inform  the  user  of  their  presence.  This  involves  either  providing 
facilities  to  correct  the  error  and  continue  if  possible,  or  executing  a 
graceful  program  exit  in  the  event  of  fatal  errors.  Meaningful  error 
messages  must  be  displayed  in  order  that  proper  corrective  action  can  be 
taken. 

2.  Program  Control 

User-friendly  techniques  can  also  be  easily  applied  to  interactive 
control  during  program  execution.  A  popular  method  to  do  this  involves 
menu-driven  selections  where  various  choices  are  displayed  along  with 
simple  commands.  The  user  selects  a  command  thereby  activating  either 
the  desired  program  segment  or  additional  sub-menus  as  appropriate.  By 
using  single  keystroke  commands,  combined  with  range  screening  and 
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suitable  error  messages,  a  very  large  number  of  alternatives  can  quickly 
be  considered. 

Displays  should  be  designed  so  that  the  user  has  enough  informa¬ 
tion  to  make  decisions  without  being  overburdened  visually  by  extraneous 
material.  Selections  should  be  arranged  with  the  most  frequently  used 
choices  the  easiest  to  invoke.  A  good  example  of  this  is  providing 
access  to  default  parameters  only  on  demand.  These  values  are  then 
transparent  to  the  user,  yet  readily  available  if  changes  are  to  be 
made. 

3.  Programming  Tactics 

The  addition  of  user-friendly  features  invariably  results  in 
more  complicated  code  and  longer  development  time.  This  is  a  small  price 
to  pay  in  view  of  the  utility  gained.  Programs  that  are  difficult  to  use 
typically  go  unused  and  often  must  be  changed,  at  far  greater  cost  than 
an  equivalent  original  design,  to  instill  user  confidence.  If  user- 
friendly  facilities  are  incorporated  at  the  conceptual  phases  of  software 
development  and  sensibly  blended  into  the  program  structure,  their  cost 
can  be  minimized. 

For  optimization  programs,  there  is  a  clear  delineation  between 
interface  and  solution  modules.  Even  with  time-shared,  or  dedicated 
microcomputer  systems,  the  code  supporting  user-friendly  facilities  can 
be  easily  isolated  from  the  iterative  portions  of  the  program.  In  this 
manner,  the  programmer  can  take  advantage  of  these  indispensable  input/ 
output  techniques  without  impairing  the  solution  process.  With  proper 
memory  management  discipline,  the  only  penalty  is  increased  secondary 
(offline)  storage  consumption. 


36 


III.  ALGORITHMS 


The  solution  programs  chosen  for  inclusion  in  this  package  represent 
advanced  methods  for  solving  minimum  cost  network  flow  problems.  Detailed 
descriptions  of  the  algorithms  employed  are  scattered  among  various 
references,  and  are  not  available  in  a  collected  form.  This  section 
outlines  the  fundamental  ideas  underlying  each  algorithm  to  provide  a 
single  source  document.  The  discussions  are  necessarily  brief  and  the 
reader  is  directed  to  the  primary  references  for  in-depth  descriptions. 

A.  GNET 

GNET  is  an  extremely  efficient  and  elegant  code  that  solves  network 
flow  problems  with  linear  costs  and  bounded  variables.  GNET  uses  the 
well-known  primal  revised  simplex  algorithm  specialized  for  networks. 

1.  Primal  Revised  Simplex  Algorithm 

Consider  the  following  linear  program: 

(LP)  MIN  cTx 

s.t.  Ax  =  b 

0  <  x  <  u 

where  c  is  a  vector  of  cost  coefficients  and  A  is  a  matrix  of  technological 
coefficients.  Any  lower  bounds  on  the  variables,  x,  have  been  eliminated 
by  transformation.  The  matrix  A  may  be  partitioned  into  two  sub-matrices 
B  and  N  to  yield 

A  ■  [B  :  N] 

where  a  matrix  of  linearly  independent  columns,  a  basis,  is  represented 
by  B  and  N  consists  of  the  remaining  columns  of  A. 
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The  vectors  c  and  x  may  be  similarly  partitioned  to  form 
^  3  CN^  ’ 

^  ”  ( Xg  Xj^ )  , 

which  implies  that 

C  X  =  C  Xg  +  c  x^  . 

Redefining  variables  so  that  every  nonbasic  variable  is  at  its  lower 
bound  simplifies  the  procedure.  Utilizing  the  notation  of  Bradley,  Brown 
and  Graves  [Ref.  3],  let  xk  be  replaced  by  uk  -  x^  (reflection)  whenever 
x^  reaches  its  upper  bound. 

Given  a  feasible  basis,  there  exists  a  unique  solution  x,  such  that 
Bx  *  b,  and  a  basic  solution  where  x^  =  (xg  5)  . 

Any  solution  satisfying  the  constraints  can  be  written  as 
Ax  =  [B  N](xB  xn)T  =  Bxg  +  NxN  =  b  , 


or 

xB  *  B-1b  -  B“^NxN  . 

Denoting  B’*N  as  Z  yields 

xB  =  B-1b  -  ZxN  . 

The  value  of  the  objective  function  c^x,  expressed  in  terms  of  xN  is 
then  cTx  =  CgB_1b  +  (cj  -  CgB’*N)xN  . 

Differentiating  with  respect  to  x^  provides  the  rate  of  change  of  the 
objective  function  in  response  to  changes  in  x^: 


acTx  _  -T  -V1M 
Tx^  CN  "  cBB  N  ’ 


The  vector  CgB’^  is  the  solution  to 
wTB  =  c l  , 


(1) 


where  w  is  commonly  referred  to  as  the  dual  solution  vector. 
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Favorable  movement  of  the  objective  function  (minimization) 
from  a  feasible  solution  is  indicated  by 

c  ■  -  <  0  (2) 

J  J 

for  nonbasic  variables  j  at  zero  flow.  Any  nonbasic  variable  x. 

J 

satisfying  (2)  will  induce  the  following  change  in  the  solution  (assuming 
all  other  nonbasic  variables  remain  fixed): 

Xg  =  B_1b  -  Z^Xj  (3) 

x^  -  ( 0 , . . • , x j , . . .  0)  .  (4) 

As  this  solution  satisfies  the  explicit  constraints,  enforcement  of  the 
bound  constraints  will  ensure  that  feasibility  is  maintained.  Updating 
the  variable  partition  completes  the  procedure. 

Bradley,  Brown,  and  Graves  [Ref.  3]  give  the  following  interpre¬ 
tation  of  the  revised  simplex  procedure: 

STEP  0:  Obtain  a  feasible  solution. 

REPEAT 

STEP  1:  Priceout.  Select  a  candidate  variable  to  enter  the 
basis  that  satisfies  (2). 

STEP  2:  Ratio.  Find  the  greatest  bound  such  that  the  incoming 
vari able: 

a.  does  not  exceed  its  upper  bound,  or 

b.  drives  a  basic  variable  to  its  lower  bound,  or 

c.  drives  a  basic  variable  to  its  upper  bound. 

STEP  3:  Pivot.  Update  the  solution  using  (3)  and  (4).  If 

case  (a)  of  STEP  2  applies,  reflect  the  candidate  incoming 
variable  chosen  in  STEP  1  and  leave  the  basis  and  dual 
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solution  unchanged.  For  case  (b),  change  the  basis 
(pivot)  and  find  a  new  dual  solution.  Case  (c)  requires 
that  the  outgoing  variable  first  be  reflected  and  then  a 
basis  and  dual  solution  update  performed. 

UNTIL  STEP  1  fails  to  find  a  favorable  candidate  (optimality). 

2.  Network  Specializations 

Recall  NPP3 

(NPP)  MINIMIZE  f (x ) 
s.t.  Ax  =  b 

T  £  x  <_  u 

where  A  is  a  node-arc  incidence  matrix.  Restricting  f  (x)  to  be  a 
linear  function  yields  a  linear  network  programming  problem  (LNP) : 

(LNP)  MINIMIZE  cTx 

s.t.  Ax  *  b 

T  £  x  <  u 

Two  well-known  results  that  characterize  the  bases  of  LNP  allow 
for  an  extremely  efficient  specialization  of  the  primal  simplex  procedure 
to  be  applied  to  LNP.  First,  all  bases  for  LNP  are  formed  by  a  set  of 
columns  which  correspond  to  a  spanning  tree  for  the  graph  represented 
by  LNP.  Additionally,  any  such  basis  matrix  B  can  always  be  placed  in 
triangular  form  by  simple  row  and  column  permutations.  The  background 
pertaining  to  these  results  can  be  found  in  any  elementary  network 
programming  reference  (e.g.,  Kennington  and  Helgason,  [Ref.  26]). 

Tri angulation  of  the  basis  matrix  B  ’"njlies  that  Z,  which  is 
defined  as  B”*N  and  thus  the  solution  to  BZ  =  N  ,  can  be  obtained 

3 See  Section  I. 
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by  direct  solution  (back  substitution).  Also,  the  dual  variables  w 
can  be  found  by  forward  substitution  of  (1). 

Characterization  of  network  bases  represented  as  spanning  trees 
for  the  underlying  graph  suggests  an  easy  method  to  perform  the  ratio 
test  and  pivot  steps  of  the  simplex  method.  Since  a  spanning  tree  by 
definition  [Ref.  26]  is  a  connected  acyclic  graph,  any  incoming  nonbasic 
variable  will  form  a  cycle  with  the  basis  tree.  Changes  in  flow  as  a 
result  of  an  introduced  arc  will  exclusively  occur  in  unit  amounts  on  the 
cycle  in  question.  All  off-cycle  basic  arcs  will  be  unaffected.  The 
direction  of  flow  change  (decrease  or  increase)  will  depend  on  the 
orientation  of  a  particular  basic  arc  with  respect  to  the  incoming  arc. 
The  basic  arc  on  the  cycle  that  reaches  its  bound  first  is  removed  from 
the  basis  (unless  the  incoming  arc  reaches  its  bound  first). 

Bradley,  Brown,  and  Graves  [Ref.  3]  employ  the  following 
algorithm: 

STEP  0:  Given  a  feasible  starting  solution 
REPEAT 

STEP  1:  Priceout.  Given  the  dual  variables  w,  the  priceout 
becomes  c^  -  w^  +  Wj  for  nonbasic  arc  k  that  is 
directed  from  node  i  to  node  j.  Select  an  incoming  arc 
k  with  c^  -  w^  +  Wj  <  0  or  terminate  with  the 
current  solution  as  optimal. 

STEP  2:  Ratio  Test.  Three  cases  analgous  to  bounded  variable 
simplex  for  which: 


a.  The  incoming  arc  reaches  its  opposite  bound  u^. 
Otherwise  a  basic  variable  is  selected  as  the  outgoing 
variable  (for  basic  arcs  on  cycle  with  incoming  arc  with 
opposite  orientation). 

b.  Incoming  arc  drives  basic  variable  down  to  its 
lower  bound,  (for  basic  arcs  on  cycle  with  incoming 
arc  with  same  orientation). 

c.  Incoming  arc  drives  basic  variable  up  to  its  opposite 
bound. 

STEP  3:  Update  (depending  upon  ratio  test  case  result). 

a.  Reflect  incoming  arc. 

b.  Pivot  update  (see  below). 

c.  Reflect  outgoing  arc,  pivot  update. 

Pivot  Update:  simultaneously  perform  a  one-pass 
update  to: 

1)  modify  flow  on  arcs  in  cycle  by  constant  equal 
to  ratio  result, 

2)  update  basis  tree  representation  ("rehang"), 

3)  update  dual  variables  for  nodes  whose  precedessor 
chain  to  the  root  is  changed, by  pivot  operation 
(these  are  the  rehung  nodes). 

UNTIL  OPTIMAL. 

The  algorithm  requires  only  integer  arithmetic  with  most  operations 
involving  only  addition  and  subtraction.  Use  of  programming  tactics  and 
coding  efficiencies  described  in  Reference  3  results  in  an  extremely 
efficient  and  compact  code,  ideal  for  microcomputer  applications. 
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B.  NLPNET 


Minimum  cost  nonlinear  network  flow  problems  are  solved  utilizing 
NLPNET  [Ref.  19],  a  primal  approach  based  on  controlled  truncation  of  a 
conjugate-gradient  method  for  solving  the  Newton  equations.  As  such, 
NLPNET  is  an  extension  of  the  unconstrained  Truncated  Newton  methods 
described  in  [Ref.  27].  Operating  on  a  maximal  basis  in  the  manner  of 
Dembo  and  Klincewicz's  Scaled  Reduced  Gradient  method  [Ref.  28],  NLPNET  is 
designed  to  take  advantage  of  the  highly  desirable  local  convergence 
properties  afforded  by  Newton  methods  while  avoiding  the  global  conver¬ 
gence  problems  and  computational  overhead  traditionally  associated  with 
these  procedures.  These  methods  belong  to  the  class  of  "inexact  Newton 
Methods"  [Ref.  29]. 

1.  Basic  Properties  of  Solutions  and  Methods  for  Nonlinear 

Optimization 

Consider  the  unconstrained  nonlinear  programming  problem  (NLP). 

(NLP)  MINIMIZE  f ( x ) 

x  c  Rn 

where 

f:  Rn*R  is  a  nonlinear  function  with  the  following  properties. 

a.  f  is  continuous  and  twice  differentiable. 

b.  At  a  relative  minimum  x*,  the  gradient  g(x*)  vanishes  and  the 
hessian  matrix  H(x*)  is  positive  definite. 

c.  For  all  x,cRn,  the  level  sets 

L(x«)  s  [x  :  f(x)£  f(x,)]  are  bounded 
(e.g.,  Oembo  and  Steihaug,  [Ref.  27]). 
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Algorithms  to  solve  NLP  vary  considerably  in  specific  approach, 
however,  most  are  iterative  descent  methods.  (Iterative  means  that  the 
algorithms  generate  a  series  of  points  based  on  the  preceeding  points 
while  descent  implies  that  each  new  point  produced  by  the  algorithm 
improves  the  solution  by  reducing  the  value  of  the  objective  function.) 

In  this  manner,  such  methods  would  in  the  ideal  case  converge  to  a 
solution  point  of  NLP. 

For  NLP,  we  must  distinguish  between  two  types  of  solution  points: 
local  minimum  points  and  global  minimum  points. 

Definition  [Ref.  30].  A  point  x*  e  Rn  is  a  local  minimum  point  of 
NLP  if  there  is  some  c  >  0  such  that  f(x)  f(x*)  for  all  x  e  Rn  within 
distance  c  of  x*. 

Definition  [Ref.  30].  A  point  x*  e  Rn  is  a  global  minimum  point  of 
NLP  if  f ( x )  >  f ( 5*)  for  all  x  e  R. 

Similarly,  the  concept  of  convergence  can  be  viewed  in  a  global 
and  local  context.  Global  convergence  deals  with  the  actual  determina¬ 
tion  of  a  solution  point  from  an  arbitrary  starting  point  and  is  not  to  be 
confused  with  convergence  to  a  global  optimum. 

Definition  [Ref.  30].  An  algorithm  is  said  to  be  globally 
convergent  if  for  arbitrary  starting  points,  it  is  guaranteed  to  generate 
a  series  of  points  converging  to  a  local  solution  point  (local  minimum). 

Local  convergence,  on  the  other  hand,  refers  to  the  properties  of 
the  algorithm  in  the  vicinity  of  a  solution.  For  example,  asymptotic 
speed  of  convergence.  One  of  the  more  important  measures  of  this  speed 
is  order  of  convergence. 
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Definition  |.Ref.  30].  Let  the  sequence  (r^l  converge  to  r*. 
The  order  of  convergence  of  {r^}  is  defined  as  the  supremum  of  the 
non-negative  numbers  p  satisfying 


0  <  lim 


I  Vi  •  r5 ' 


k“  I  r. 


-  r* 


(5) 


Larger  values  of  p  in  (5)  imply  faster  convergence  since  the  distance 

1  L 

from  the  limit  r*  is  reduced  by  the  pl  power  in  a  single  step. 

For  example,  quadratic  convergence  (p  =  2)  doubles  precision  at  each 
step  while  for  linear  convergence  (p  =  1)  the  reduction  ratio  remains 
constant. 

Finally,  if  a  sequence  has  order  p  convergence,  then  as  k  becomes 
very  large,  we  have 


1  Vl  -  r*  I  =  9  I  rk  .  r*  \*  , 
where  8  is  a  constant  termed  the  convergence  ratio  [Ref.  30]. 


2.  Newton's  Method  for  Unconstrained  Optimization 

NLP  can  be  solved  using  the  well-known  Newton  method 
[e.g..  Ref.  30],  where  the  function  of  f  being  minimized  is  approximated 
locally  by  a  quadratic  function  and  this  approximating  function  is  then 
minimized  exactly.  Thus  near  the  point  x0,  we  can  approximate  f  by  the 
truncated  Taylor  series: 

f (x)  *  f(x,)  +  g(x,)T(x  -  xf)  +  l/2(x  -  x,)TH(x,)(x  -  x,)  . 


iX’+'r- *>* !*'***<  •• 
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Minimizing  the  approximation  by  setting  the  derivative  equal  to  zero 
yields: 

g'(x)  =  g(x,)  +  g'(x,)(x  -  x0)  =  0  , 

where  g(x#)  is  defined  as  the  gradient  vector  of  the  function  f 
evaluated  at  the  point  x,  and  g'(x,)  is  equivalent  to  H(x,),  the 
hessian  matrix  of  f  evaluated  at  x#. 

Making  the  appropriate  symbol  replacements  and  defining  the 
vector  (x  -  x,)  as  p  gives  us: 

H(x,)  p  =  -g(x.)  .  (6) 

Equation  (6)  is  the  classical  representation  of  the  Newton  method  in 
an  unconstrained  setting  with  p  often  referred  to  as  the  Newton 
direction  or  the  Newton  step.  This  suggests  the  following  iterative 
procedure  given  an  initial  guess  x0: 

Step  0.  k  =  0. 

REPEAT 

Step  1.  Solve  H(xk)pk  =  -g(xk)  Tor  pk  . 

Step  2.  xk+1  *  xk  +  pk  . 

UNTIL  convergence. 

An  equally  familiar  result  is  that  under  the  regularity  assump¬ 
tions  given  earlier  this  method,  although  well  defined  near  a  solution 
point,  may  not  converge  when  far  from  a  solution.  One  reason  for  such 
behavior  is  that  the  quadratic  function  is  a  poor  approximation  of  f  at 
the  point  xk>  Since  in  its  pure  form  the  Newton  direction  contains 
both  direction  and  steplength  (the  magnitude  of  the  direction  vector). 
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information  based  on  the  assumption  of  quadratic  properties  for  f,  either 
or  both  of  these  quantities  may  be  inappropriate  for  the  function  under 
consideration.  Traditionally,  the  method  has  been  modified  to  enhance 
its  convergence  properties  and  accommodate  arbitrary  functions.  The  most 
frequent  modification  is  the  introduction  of  a  search  parameter  x  as 
follows: 


Vi  *  *k *  Vic  •  (7) 

where  x^  is  chosen  to  minimize  f  in  the  direction  p^.  For  points 
where  the  Newton  quadratic  approximation  is  "good,"  xk  should  be 
approximately  one.  When  xfc  i  1,  the  ray  search  (7)  is  used  as 
insurance  that  descent  with  respect  to  f  is  realized  on  the  Newton 
direction. 

Finally,  Newton's  method  can  fail  if  H(xk)  becomes  non-positive 
definite  at  points  far  from  a  solution.  In  such  cases  the  method  may  not 
yield  a  direction  of  descent  and  can  in  fact  seek  a  relative  maximum 
point.  The  method's  characteristics  are  sunmarized  in  Figure  4. 

3.  Truncated-Newton  Methods  for  Unconstrained  Optimization 

Since  the  benefits  of  the  Newton  direction  are  greatest  in  the 
immediate  vicinity  of  a  solution,  where  the  quadratic  approximating 
function  best  describes  the  actual  function,  an  accurate  determination  of 
the  Newton  Direction  appears  unnecessary  when  far  from  a  solution.  The 
use  of  iterative  methods  to  solve  the  Newton  equations  suggests  a  direct 
trade-off  between  computational  effort  and  Newton  solution  accuracy.  It 


ADVANTAGES 


1.  Locally  and  quadrat i cal  1 y  convergent. 

DISADVANTAGES 

1.  The  method  is  not  globally  convergent. 

2.  It  is  not  defined  at  points  x,  where  H(x) 
is  singular  or  essentially  singular  in  a 
numerical  sense. 

3.  For  nonconvex  problems,  it  does  not 
necessarily  generate  a  sequence  of  descent 
directions. 

4.  An  n-dimensional  linear  system  of. equations 
must  be  solved  at  each  iteration. 


Fig.  4.  Characteristics  of  the  Newton  Method. 

is  precisely  this  relationship  that  the  Truncated-Newton  class  of  methods 
[Ref.  27]  seeks  to  exploit.  At  the  outset  of  the  optimization,  easily 
obtainable  but  relatively  inaccurate  approximations  to  the  Newton  Direc¬ 
tion  are  tolerated  with  increasing  accuracy  demanded  as  x^  approaches  x*. 

In  order  to  control  such  a  method,  a  measure  of  required  accuracy 
is  needed  that  reflects  how  "far"  the  present  value  of  the  objective 


function  is  from  a  solution.  The  currently  preferred  measure  [Ref.  27] 
is  the  relative  residual 

II  r^  ||  /  ||  ||  where  r^  is  defined  as 

\  =  H< VP|c  +  9(*k)  • 

The  iterative  method  applied  to  the  Newton  equations  is  therefore 
truncated  when  the  relative  residual  is  "small  enough." 

The  basic  outline  of  a  Truncated-Newton  method  [Ref.  27]  is 
given  in  Figure  5.  Equations  (8)  and  (9)  are  conditions  guaranteeing 
"sufficient  descent"  [Ref.  27].  It  can  be  shown  [Ref.  31]  that  if  these 
conditions  are  satisfied  for  x^  and  if  {xk>  converges  to  an  optimal 
point  x*  at  which  H(x*)  is  positive  definite,  then  there  is  an  iteration 
index  k  >  0  such  that  x^  =  1  admissible  for  k  _>  k,  and  {xk>  converges 
superlinearly  (order  >  1)  to  x*.  These  conditions  therefore  operate  in 
conjunction  with  the  direction  finding  mechanism  to  produce  desirable 
terminal  convergence  properties  (when  the  method  does  in  fact  converge). 

Using  a  conjugate-gradient  (CG)  algorithm  to  iteratively 
calculate  the  search  direction  pk,  Dembo  and  Steihaug  [Ref.  27]  pro¬ 
pose  the  truncated  conjugate-gradient  algorithm  (TNCG)  illustrated 
in  Figure  6  to  serve  as  the  minor  iteration. 

There  are  three  different  ways  in  which  the  TNCG  minor  iteration 
can  terminate: 
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Given  an  initial  solution  x#: 


k  ♦  0 
REPEAT 

!  Major  J - 

I  Iteration  j 

I  Compute  f(xk),  g(xk) ,  H(xk) 

j  Test  for  convergence 

|  IF  (NOT  convergence)  THEN 

I  _ 

j  J  Minor  J 

|  j  Iteration  j - 

|  |  Calculate  pk  such  that 

I  I  H(xk)pk  ■  -g(xk) 

I  j  Using  a  Truncated  Newton  iterative  method 

I  I _ 

I 

I  Calculate  some  xk  that  satisfies 

j  f(xk  +  ^kPk)  1  f(xk)  +  aXic9(xk)Tpk  ;  a  e(0,l/2)  (8) 

I  .  .  T-  _  .  T- 

I  g(xk  +  xkpk)  Pk  _>  eg(xk)  Pk  ;  8  e  (a,  1)  (9) 

J  Vi  *  x'k  +  xkPk 

I  k  *  k+1 


END  IF 


Until  convergence 


Fig.  5.  Truncated-Newton  Method 
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Denoting  k  as  the  major  iteration  counter  and  f(xk),  g^),  H(x^) 


as  fk,  gk  and  respectively 


Minor 


,  ninur  - 

Iteration 


Step  1:  p,  *  0 

-  -gk 
a,  «■  r. 


REPEAT 
Step  2: 


Step  3: 


aTdt  ♦  r 1?  , 


q,  .  Hka( 

IF  aTq1  £  c  (ajai )  THEM 


EXIT: 


END  IF 


3#  if  i  =  0 


p.  otherwise 


a.  ♦  r^r .  /  a^q. 

Sin  *  Pi  +  °iai 


rin  *  ri  -  Vi 

IF  II  rj+1  II  /  II  gk  II 

Em:  pk  *  pi+1 

ENOIF 

Step  4:  s.  -  rT+1  fj+1  /  i=Tf . 

ain  *  r'm  *  Vi 
ainam  *  finfin  *  9iaIai 


<  nk  THEN 


i  *  i+1 


UNTIL  EXIT 


Fig.  6.  Truncated  Conjugate-Gradient  Algorithm 


a.  The  gradient  vector  gk  points  in  a  direction  of 
negative  curvature  (g^H^g^  <  0).  In  this  case,  the  minor  iteration 
returns  with  3  -g^,  the  steepest  descent  direction. 

b.  A  direction  of  negative  curvature  is  encountered  in  the 

T  * 

CG  iteration  (3-Hk3-  <  0)  prior  to  satisfying  the  Truncated-Newton 
termination  criterion.  The  CG  procedure  is  terminated  and  the  current 


estimate  pk  is  used.  Oembo  and  Steihaug  [Ref.  27]  show  that  p^  is  a 
direction  of  descent. 


c.  The  algorithm  terminates  with  the  Truncated-Newton 
criterion.  Dembo  and  Steihaug  [Ref.  27]  show  that  this  always  occurs 
in  the  vicinity  of  a  strong  local  minimum  and  thus  is  the  determining 
factor  in  the  rate  of  convergence. 

Dembo  and  Steihaug  [Ref.  27]  define  {n^l  as  the  forcing 
sequence,  case  c  above  as  the  Truncated-Newton  termination  and  directions 
resulting  from  either  cases  a  or  c  as  Truncated-Newton  directions. 

It  can  be  shown  [Ref.  27]  that  the  TNCG  algorithm  is  globally 
convergent  and  capable  of  coping  with  regions  where  the  hessian  is 
indefinite.  The  following  theorem  indicates  that  the  order  of  convergence 
can  be  controlled  with  the  proper  choice  of  the  forcing  sequence. 

Theorem  1.  (Dembo  and  Steihaug  [Ref.  27])  Properties  of  the 
Forcing  Sequence. 

Assume  that  the  Truncated-Newton  iterates  {x^}  converge  to 
x* .  Then 
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a.  ix^}  ♦  x*  superlinearly  (order  >  1)  if  {n^}  +  0. 

b.  (xk)  +  x*  with  order  (1+t)  if 

nk  1  c  11  9k  1 

for  some  c  >  0  and  0  <  t  <  1. 

c.  {xk}  x*  linearly  (order  1)  if  nk  is  uniformly  less 
than  one  and  bounded  away  from  zero. 

Theorem  1  indicates  that  by  choosing 

nk  =  c  M  gk  \\l  (k  =  0,1,2,...) 

for  some  c  >  0  and  0  <  t  <_  1,  the  TNCG  algorithm  will  possess 
any  prescribed  order  of  convergence  (1+t)  between  one  and  two. 

The  result  is  an  adaptive  algorithm  that  solves  NLP.  When  far 
far  from  a  solution,  ||  gk  ||  is  large  as  is  nk  and  little 
effort  is  needed  to  satisfy  the  Truncated-Newton  criterion.  As 
(V  +  x*,  {gk>  ♦  0  which  implies  {nk>  +  0  thereby  forcing  {?k>  ♦  0  and 
Pk  to  the  Newton  direction.  The  method  therefore  automatically 
incorporates  an  increasing  amount  of  second-order  information  as  the 
optimization  process  progresses  just  when  such  information  is  of 
greatest  use. 

4.  Truncated-Newton  Methods  for  Linearly  Constrained  Optimization 

Truncated-Newton  methods  derived  for  unconstrained  optimization 
can  be  extended  to  solve  linearly  constrained  nonlinear  problems  of  the 
form 
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(LCNLP)  MINIMIZE  f(x) 
x  e  Rn 

s.t.  A  x  =  b 
1  £  x  _<  u, 

where  the  function  f { x )  is  convex  and  separable  (i.e., 
f ( x )  =  f .(x.)  and  each  f,(x.)  is  a  convex  function). 

J  "*  J  J  J  J 

a.  The  "Reduced  Problem" 

In  the  manner  of  Murtagh  and  Saunders  [Ref.  32],  partition 
the  columns  of  A  as  follows: 

A  =  [B  S  N], 

where  the  columns  of  B  form  a  basis;  the  columns  of  S  correspond  to  super- 
basic  variables  (nonbasic  variables  whose  flow  is  allowed  to  vary 
between  bounds);  and  the  columns  of  N  to  nonbasic  variables  with  flow 
fixed  at  either  bound  and  not  allowed  to  vary.  Similarly,  the  other 
important  vectors  can  be  partitioned. 

x  =  [xB  xs  J/  (primal  solution), 

g(x)  ■  Cg(x8)  g(x$)  g(xN)]  (gradient  vector), 

p  «  [pg  p<-  pN]  (search  direction  vector), 

as  can  the  function  f(x)  =  f(Xg,Xg,xN). 

This  partition  allows  us  to  re-express  the  constraint  set 

A  x  =  b 
as 

[B  S  N][xb  x$  xN]  3  b 
or 

B  Xg  +  S  x$  +  N  xN  =  b 
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solving  for  xg  (utilizing  the  fact  that  B  represents  a  basis  and  thus 
B’*  exists), 

xB  =  B-1b  -  B-1Sxg  -  B_1Nxn  . 

Now  we  are  able  to  express  f  completely  in  terms  of  [x^  x^]  and  denote 

A 

this  function  by  f. 

f(xs,  xN)  =  f(B_1b  -  B‘1Sxs  -  B_1Nxn,  xs.  xn)  . 

a 

Calling  f(x^,  xN)  the  "reduced  problem"  or  RNLP,  we  note  the  following 


g  -  f'(5s.  V  *  [gs  9N] 

and  by  the  chain  rule, 

A 

r  _  3f  .  3f  3f  3Xn  _  - 

9S  '  lxs  ‘  lxs  TXg  15^  ‘  9S  * 

A 

z  _  af  _  af  .  af  axj,  _  - 

9N  "  TxN  ‘  lxN  TxB  “  gN  * 

Therefore  the  modified  problem  becomes 

(RNLP)  MINIMIZE  f(x$,  xN) 

x$,  xN  e  R 

s.t.  S.  *s  — 

—  XN  —  UN  * 


gB(B‘1S) 
gB(B"1N)  . 


(10) 

(ID 


As  the  bound  constraints  (10)  and  (11)  can  be  handled  implicitly,  we 
essentially  are  now  dealing  with  an  unconstrained  problem, 
b.  The  Search  Direction 

A  generic  minimization  algorithm  for  RNLP  would  calculate 
successive  search  directions  from  a  feasible  initial  point,  updating  x 
until  convergence.  Since  for  iteration  k  of  this  sequence, 

“k  '  Ok  *  Vi  -  ;k 
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and  from  previous  results  we  know  that 
axk  *  [axB  ax$  axN]  , 


axN  3  6  by  definition. 


axB  3  B-1b  -  B^Sxg 


(k+1) 


B_1Nxw  -  B_1b  +  B_1Sx,  +B"1Nxm 
N(k+1)  b(k)  V) 


3  -B'1S(x«.  -  x,  )  3  -B-1Sax,  . 

a(k+l)  a(k)  5 


It  follows  that 


ax 


(k) 


axr 


ax. 


l(k) 

’(k) 


L%)J 


“  - 1 
-B  Sax, 

-B-1S 

b(k) 

““ 

ax, 

b(k) 

I 

0 

.  0  . 

ax 


sik)  • 


and 


Pk  *  **k 


-B_1S 

I 

0 


Ps 


(k) 


Given  a  feasible  point  x#,  a  Newton  method  to  solve  RNLP 
would  involve  successive  solutions  of  a  quadratic  program  to  compute 
search  directions. 

MINIMIZE  l/2pTH(x#)p  +  g(x.)Tp  (12) 

P 


s.t.  Pj<0, 
Pj  >  °, 


Denoting  the  matrix 


lj  U  *  S) 

Uj  (j  e  B)  . 

as  Z  and  altering  (12)  for  RNLP  implies 
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MINIMIZE 

h 

l/2pJzTH(;.)Zps  +  g(x,)TZps 

(13) 

s.t. 

PJ  i  °-  xj  =  'j 

(14) 

Pj  >  0,  Xj  »  Uj  . 

(15) 

The  feasibility  constraints  (14)  and  (15)  limit  admissible 
search  directions  for  variables  at  bounds.  As  pN  is  chosen  to  be 
zero  and  can  be  handled  explicitly  during  solution  of  (13),  only 
the  components  of  pg  are  likly  to  cause  trouble.  Dembo  and  Klincewicz 
[Ref.  28]  point  out  that  the  major  difficulty  is  one  of  wasted  computational 
effort.  If  after  p^  is  determined  the  induced  pB  violates  (14)  or  (15), 
a  new  basis  must  be  found  and  a  new  search  direction  [jjg 
until  all  constraints  are  satisfied.  They  further  show  that  if  one 
operates  using  a  basis  with  the  greatest  number  of  variables  between 
bounds  (free  variables),  the  only  nonbasic  variables  that  cause  concern, 
with  respect  to  feasibility,  are  those  moving  away  from  a  bound.  Such  a 
basis  is  termed  a  maximal  basis  [Ref.  28]. 

The  solution  to  (13)  when  H(x#)  is  positive  definite  (recall  f 

A 

and  thus  f  are  convex)  is, 

ZTH(x,)Zps  =  -ZTg(xf )  . 

ZTH(x,)Z  and  ZTg(x,)  are  often  referred  to  as  the  reduced  hessian  and 
and  reduced  gradient  respectively. 

We  now  have  a  system  of  equations  in  the  unknown  vector  p 
that  provides  the  solution  to  the  direction  finding  problem  for  RNIP 


p<.]  calculated 
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and  hence  for  its  equivalent,  NLP.  Additionally,  we  can  recover  the 
complete  p  vector  through  the  relationship, 

p  =  Zp$  . 

The  Truncated-Newton  Congugate-Gradient  method  can  now  be  applied 
directly  to  the  quadratic  approximation  of  RNLP. 

5.  Primal  Truncated-Newton  Method  Specialized  for  Networks 

The  material  presented  to  this  point  is  totally  general  and 
nothing  in  its  development  relies  on  the  fact  that  the  problem  in  question 
is  a  network  program.  With  this  in  mind,  it  is  now  time  to  examine  some 
of  the  special  structure  afforded  by  the  network  representation. 

Define  A  of  LCNLP  (and  RNLP)  to  be  a  node-arc  incidence  matrix 
of  a  directed  network.  The  result  is  the  usual  formulation  of  the 
capacitated  minimum  cost  network  problem  with  a  convex  separable  objective 
function.  We  shall  label  these  new  problems  NLN  and  RNLN. 

Any  solution  technique  for  NLN  can  naturally  take  advantage  of 
the  special  structure  of  network  optimization  problems,  dispense  with 
explicit  representations  of  the  basis  inverse  and  perform  update  and 
solution  operations  directly  on  a  network  specialization  of  the  basis 
representation. 

First  consider  the  form  of  the  reduced  gradient  of  RNLN, 

ITg  *  [-sVT  1  0]  [gB  g$  g„]T 
■  -sVTg„  *  5S  . 

Recognizing  B"^gg  as  the  estimates  of  the  dual  variables  (denote 
as  w)  and  the  solution  to 


58 


suggests  that  w  can  be  determined  by  solving  a  triangular  system  of 
equations  since  network  bases  B  (and  hence  their  transpose)  can  always  be 
placed  in  triangular  form  by  simple  permutation.  Calculation  of  the  dual 
variables  in  this  manner  would  allow  easy  computation  of  the  reduced 
gradient.  This  is  a  relatively  standard  approach  to  NLN. 

The  solution  technique  just  developed  serves  as  a  framework 
in  which  to  embed  an  adaptive  direction  finding  mechanism  and  thus 
produce  a  Truncated-Newton  method.  With  the  addition  of  control  logic 
and  procedures  to  monitor  and  manipulate  the  variable  partition,  we 
have  all  the  necessary  ingredients  of  a  complete  algorithm  to  solve 
NLN  [Ref.  19]. 

Consider  NLN  and  the  associated  problem  RNLN.  Given  an  initial 
feasible  solution  x,, 

STEP  0:  Partition  x,  into  basic,  superbasic  and  nonbasic  sets  in 

such  a  manner  that  a  maximal  basis  is  established.  Likewise, 
partition  g  and  H. 
k  «■  0 

STEP  1:  Calculate  the  dual  variable  estimates  (w)  by  solving 
BTw  *  -gB  (w  *  -B"TgB). 

STEP  2:  Compute  the  reduced  gradient, 

lrs  *  zT9  *  -sTb'T9b  ♦  9$  ■  STi  *  gs  . 
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STEP  3:  Test  for  optimality  on  subspace  defined  by  active  superbasic 


variables. 

IF  | |  Z^g  | |  <  tolerance  THEN  optimal  on  current  subspace 

FIND  those  nonbasic  variables  eligible  to  enter  superbaisc 
set.  Eligibility  is  determined  by  potential  for  feasible 
displacement  from  bound  and  subsequent  reduction  in  objective 
function  value. 

A 

Compute  *  -NTB"TgB  +  9N  =  NTw  +  §N  . 

N 

a 

IF  (sf/sx^  0  and  (xN)j  =  lj)  OR 

A 

(3f/axN_<  0  and  (xN)j  -  Uj)  for  all  j  €  nonbasic  set 

THEN 

STOP  -  OPTIMAL  solution  has  been  found. 

ELSE 

Add  the  nonbasic  variable(s)  not  satisfying  the  above 
conditions  to  the  superbasic  set. 

GOTO  STEP  1. 

ENDIF 


STEP  4:  Compute  a  feasible,  improving  direction  (Truncated- Newton 
direction)  by  solving 

(ZrHZ)p$  =  -ZTg 

using  a  Conjugate-Gradient  (CG)  algorithm  with  the 
termination  rule 

^  ^i  ^  £  force  . 

I  I  1%  || 

where 

forcek  =  min  (force#  ,  ||  ZTg  1 1 t > ,  force#  <  1,  t  e  (0,1], 
?i  =  -(ZTHZ(ps).  +  ZTg)  , 
i  =  CG  iteration  number. 

STEP  5:  Ensure  all  components  of  p^  calculated  in  STEP  4 
maintain  feasibility. 

Let  (ps)j  =  0  if  (p$)j  £  0  and  (x$)j  «•  1^  , 

(ps)j  *  0  if  (ps)j  £  0  and  (xs)^  *  . 

STEP  6:  Calculate  a  search  direction  for  the  basic  variables  by 
solving 

BpB  =  -Sps  . 

STEP  7:  Find  a  steplength  giving  "sufficient  descent" 
(Goldstein-Armi jo  conditions). 

STEP  8:  Update  the  flow, 

Vl  *  xk  +  Xk^k 

IF  (xb)j  at  bound  then  pivot  and  replace  with  free  arc 
from  Xg  (if  possible)  to  maintain  maximal  basis. 

IF  ( x ^ ) j  at  bound  then  remove  from  superbasic  set. 
k  ♦  k  +  1. 


GOTO  STEP  1. 


Changes  of  bases  are  performed  using  the  pivot  mechanism  of 
GNET  [Ref.  3].  The  single  variable  linesearch  to  determine  *  is  a 
safeguarded  successive  cubic  approximation  method  [Ref.  33]  modified  to 
incorporate  the  Goldstein-Armi jo  conditions,  (8)  and  (9).  This  par¬ 
ticular  linesearch  is  not  crucial  to  the  success  of  the  method  and  any 
reasonable  substitute  could  be  used;  although  quite  complex,  it  is 
reported  to  be  very  robust  [Ref.  33].  The  TNCG  algorithm  requires  that 
both  gradient  and  hessian  information  be  available.  Finally,  the  algo¬ 
rithm  requires  an  initial  feasible  solution.  This  is  not  a  limitation, 
but  rather  an  advantage  as  the  first  solution  can  easily  be  provided  by 
efficient  linear  network  programs  (i.e.,  GNET)  thereby  allowing  the 
nonlinear  code  to  be  streamlined  and  overall  computational  effort  reduced. 

C.  ENET 

ENET  solves  network  programming  problems  with  “elastic"  ranged 
constraints  and  bounded  variables  using  a  modified  revised  primal  simplex 
method . 

Elastic  constraints  can  be  violated  by  incurring  a  linear  penalty 
as  opposed  to  the  rigid  or  inviolate  constraints  of  the  classical  (network) 
programming  problem.  Brown  and  Graves  [Ref.  20]  point  out  that  such 
an  elastic  model  is  a  realistic  and  powerful  portrayal  of  many  real 
world  situations.  Hence  it  may  be  advantageous  to  relax  some  constraints 
(and  incur  a  penalty)  in  order  to  satisfy  others,  or  improve  the  value  of 
the  objective  function.  This  class  of  models,  therefore,  offers  the 
analyst  complete  flexibility  in  the  formulation  of  NPP. 


1.  Ranged  and  Bounded  Model 

The  traditional  bounded  linear  (network)  model  (LNP)4  is 
(LNP)  min  cTx 

s.t.  Ax  =  6  (16) 

1  £  x  £  u 

Addition  of  upper  and  lower  ranges  on  the  (flow  conservation) 
constraints  (16),  yield  the  following  ranged  and  bounded  model  (RLNP) : 

(RLNP)  min  cTx 

£  £  Ax  £  ¥ 
i  <  x  <  u  . 

By  introducing  additional  (structural)  variables  y  and  (logical) 
variables  s  as  follows: 

y  =  x  -  1  , 
and 

0  £  s  £  ¥  -  £  , 

RLNP  is  transformed  into  the  equivalent  equality-constrained  model  with 
translated  ranges  and  bounds: 

min  cTy  +  cTl 
s.t.  Ay  +  s  =  b  -  A1 
0  £  y  £  u  -  1 
0  £  s  £  F  ■  5  . 

s  is  a  vector  of  nonnegative  slack  variables  one  for  each  constraint, 
that  measure  deviation  of  the  current  constraint  value  from  the  appropriately 
translated  upper  range. 

^Section  III,  Subsection  A. 


2.  Elastic  Model 

Now  consider  the  enhancement  to  RLNP  where  the  constraint  (flow) 
ranges  are  allowed  to  be  violated.  For  the  purposes  of  this  paper,  the 
penalty  functions  will  be  restricted  to  the  linear  case  to  maintain 
piecewise  linearity  of  the  objective  function.  Let  7  and  ^  be  vectors  of 
the  penalty  coefficients  (the  ith  elements  define  the  cost  per  unit 
violation  for  the  ith  constraint)  for  the  upper  and  lower  constraint 
ranges  respectively.  The  resulting  model  is: 

(ERLNP)  min  c^y  +  c^l  +  7Tr  +  _zTa 
s.t.  Ay-r  +  a+  s  =  b’-Al 
6  <  y  <  u  -  1 

0<  s<b-b  (17) 

a,  r  ^  0 

The  vectors  a  and  r  are  (logical)  artificial  and  surplus  variables 
introduced  to  maintain  equality  constraints  with  all  nonnegative  variables 
as  is  the  custom  for  the  simplex  method. 

A  mathematically  equivalent  model  for  ERLNP  is 

min  c^y  +  c^i  +  z(Ay  -  b  +  Al)  (18) 

with  _  _ 

z  if  Ay  _>  b  -  Al 

z  =  1  -z  if  Ay  <  _b  -  Al 

0  otherwise  . 

7,  2  >  0  (f°r  convexity)  . 

This  is  known  as  the  Lagrangian  form,  and  illuminates  the  fact 
that  7  and  z  are  actually  bounds  on  the  variables  of  the  dual  formulation. 
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A  graphical  representation  of  the  elastic  cost  implication  for  each 
constraint  of  ERLNP  is  given  by  Figure  7. 


Cost 


bi 


Fig.  7.  An  Elastic  Ranged  Constraint 

In  this  sense,  with  i  -  7  *  «,  each  constraint  can  be  depicted  as 
in  Figure  8  where  the  discontinuities  at  the  upper  and  lower  ranges 
indicate  that  the  ranges  are  rigid  as  in  (RLNP)  and  cannot  be  violated 
(i.e.,  an  infinite  penalty  beyond  the  defined  range).  Flow  within  the 
allowable  range  incurs  zero  penalty. 
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Fig.  8.  A  Rigid  Ranged  Constraint 


When  b.  =  J)  •  =  b^ ,  the  range  collapses  to  a  point  as  shown  in 
Figure  9  and  the  model  specializes  to  LNP. 


Cost 


A,y  + 


V  5 


V 


Fig.  9.  Classical  Equality  Constraint 

Consider  the  relationship  between  the  logical  variables  s,  a,  and 
r  in  ERLNP.  The  slack  variable  si  measures  the  distance  a  solution  is 
from  the  upper  range  of  constraint,  i.  As  given  by  (17),  s^  has  an 
upper  bound  equal  to  the  difference  between  the  upper  and  lower  ranges. 

The  artificial  variable  a^  measures  the  distance  a  particular  solution 
is  below  the  lower  range  while  the  surplus  variable  r  measures  distance 
above  the  upper  range  (both  are  unbounded  variables  in  the  current 
presentation).  Figure  10  shows  this  relationship  for  a  given  constraint, 
s^,  a^,  and  r^  are  mutually  exclusive  in  any  basic  solution  since  they 
are  evidently  linearly  dependent  columns.  Additionally,  ai  and  ri  must  be 
equal  to  zero  if  non-basic,  while  s^  can  be  non-basic  at  either  zero  or  its 
upper  bound  1^  -  b^ .  Thus,  for  any  given  solution,  the  status  of  s.j, 
a1t  and  r  ■  may  be  summarized  (Figure  10): 
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1.  Upper  range  violated  by  r-  at  cost  7 ^  r ■  (a^  *  S|  =0). 

2.  No  violation  (  a^  =  r-  =  0). 

3.  Lower  range  violated  by  ai  at  cost  ai  ( s ^  *  b^  -  r-  -  0). 


Cost 


Fig.  10.  Relationship  8etween  Logical  Variables  for  ERLNP. 

3.  Elastic  Primal  Revised  Simplex  Specialized  for  Networks 
All  solutions  to  ERLNP  which  satisfy  bounds  (1 •  £  x,  _<  u.) 
are  feasible  since  the  artificial  and  surplus  variables  are  unbounded. 
Therefore,  a  solution  technique  for  ERLNP  does  not  need  to  be  partitioned 
into  the  usual  two  phases:  (I)  where  feasibility  is  achieved,  and  (II) 
where  optimality  is  obtained  while  maintaining  feasibility.  Any  starting 
solution  satisfying  bounds  1  and  u  will  serve  to  begin  the  solution 
process.  Essentially,  ENET  employs  the  primal  revised  simplex  technique  of 
GNET  with  modifications  in  logic  and  data  structure  to  accommodate 
the  additional  arcs  (variables)  implied  by  the  elastic  model.  Figure  11 
shows  the  logical  arcs  employed  for  each  node.  For  any  solution,  all 
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Artificial 


Fig.  11.  Logical  Arcs  Generated  for  Node  i 

required  information  for  each  (flow  conservation)  constraint  and  its 
penalty  contribution  to  the  objective  function  is  generated  logically 
from: 

1.  which  logical  arc  (if  any)  is  in  the  basis,  and 

2.  which  bound  the  slack  arc  s-  assumes  when  non-basic. 

Brown  and  Graves  [Ref.  20]  describe  the  following  Elastic  Primal 
Simplex  Network  Algorithm: 

STEP  0:  Select  a  starting  so!  ;fion. 


REPEAT 


STEP  1:  Price  Out.  Given  the  dual  variables  w,  the  reduced 
cost  of  each  explicit  arc  is: 


rk  -  ck 

-  W  •  +  W  . 

1  J 

(arc  k  orientated  i  to  j,  or  non- 
basic  arc  k  at  upper  bound 
orientated  j  to  i) 

and  for 

logical  arcs. 

Ui +  wi 

for  a^ 

ri  * 

0  +  wi 

for  s^ 

[Ti  ‘  wi 

for  ri  . 

Select  an  incoming 

arc  with  r.  <  0  or 

terminate  with  the 

current  solution  OPTIMAL. 

STEP  2:  Ratio  Test.  Determine  outgoing  arc  as  in  GNET. 

STEP  3:  Update.  Same  as  GNET,  substituting  obvious  specializa¬ 
tions  for  updates  involving  logical  arcs. 

UNTIL  OPTIMAL 

The  only  data  required  by  this  enhanced  elastic  model  in  addition 
to  classical  GNET  are  for  each  node  the  new  penalties,  the  difference 
between  the  upper  and  lower  ranges,  and  whether  s^  is  reflected. 

The  resulting  algorithm  and  data  structures  are  well-suited 
for  microcomputer  implementation.  The  moderate  increase  in  data  storage 
is  compensated  by  a  significant  model  enrichment.  In  particular, 

1.  Total  supplies  and  total  demands  no  longer  need  be  equal,  nor 
must  supply  nodes  be  connected  to  demand  nodes  by  paths  of 
sufficient  capacity  to  insure  classical  feasibility.  Infeasible 
problems  are  intrinsically  diagnosed  and  optimally  treated. 
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2.  All  problems  possess  feasible  primal  solutions  in  the  extended 
elastic  formulations,  and  this  provides  reliable  bounds  for  dual 
solutions. 

3.  All  bounded  problems  possess  optimal  solutions  yielding  optimal 
dual  solutions. 

4.  Informal  relaxation  (e.g.,  Lagrangian  Relaxation)  methods  are 
naturally  accommodated  in  this  context. 

5.  Formal  decomposition  methods  are  provided  much  more  robust 
primal  and/or  dual  proposals. 

6.  (Mixed)  integer  models  produce  strictly  feasible  solutions  in  this 
extended  Lagrangian  context. 

and  perhaps  most  important: 

7.  Solutions  are  reliable,  inexpensive,  and  managerially  appealing. 
Finally,  as  we  shall  see,  the  complete  Lagrangian  objective  function 
yields  a  unifying  perspective  of  solution  properties  and  intrinsically 
gives  profound  information  to  the  analyst  or  algorithms  using  ENET. 

D.  UNET 

UNET  solves  elastic  ranged  bounded  linear  network  problems  where 

certain  arcs  are  designated  as  "1-u"  or  binary  arcs.  These  arcs  can 

admit  flow  at  one  of  two  possible  values--the  lower  or  upper  bound, 

hence  the  label  "1-u".  Consider  first  the  simplest  bounded  "1-u"  model: 
(LULNP)  min  c^x' 

s.t.  Ax  =  6 

1  <  x  <  G 

x  s  {xf,  x  ),  x  c  { T ,  u >  (integer  restriction). 
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By  adding  the  now  familiar  elastic  range  framework  LULNP  becomes 
(ELULNP)  min  cTy  +  c T1  +  zTr  +  |Ta 
s.t.  Au-r  +  a  +  s  =  F-Al 
5<y<u-l 
o  <_  s  <_  F  -  5 
a,  r  _>  0 

y  =  {y^,  y^},  y  e  {T,  u)  (integer  restriction). 

Such  a  formulation  is  termed  a  mixed  integer  problem  (MIP)  since 
there  are  both  fractional  (with  respect  to  flow  bounds,  xf)  and  integral 
(again  with  respect  to  flow  bounds,  x^)  arcs  present.  Note  that  the 
cost  coefficients  cu  associated  with  yu  may  be  interpreted  as  fixed 
charges  in  the  sense  that  for  admissible  solutions,  the  cost  contribution 

-T-  -T. 

is  either  cul  or  cyu  for  each  "1-u"  arc,  and  that  each  "1-u"  arc  is 
essentially  equivalent  to  a  binary  decision  variable. 

Brown  and  Graves  [Ref.  20]  employ  an  enumeration  technique  to 
solve  ELUNLP  which  exploits  the  elastic  model  structure  and  produces 
excellent  solutions  satisfying  the  integer  restrictions  with  very  little 
computational  effort. 

Their  approach  is  analogous  to  classical  branch  and  bound  [e.g.. 
Ref.  1]  with  the  following  specializations  and  supporting  observations: 

1.  Any  continuous  relaxation  of  ELULNP  (with  integrality  violated 
by  one  or  more  arcs  in  yu)  provides  a  lower  bound  for  the 
value  of  an  optimal  solution  to  ELULNP. 

2.  Any  continuous  relaxation  of  ELULNP  can  be  rounded  to  an  integer 
solution  (satisfying  integrality  for  y  )  with  very  little 
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effort.  Further,  such  rounded  solutions  are  all  admissable 
(feasible  in  the  extended  elastic  sense). 

3.  Restrictions  of  ELUNLP  (with  one  or  more  arcs  in  y^  fixed 

at  a  bound)  are  all  admissible  (elastically  feasible)  and  possess 
solution  values  no  higher  than  the  lower  bound  of  their  respec¬ 
tive  relaxations. 

Thus,  enumeration  by  branch  and  bound  may  be  organized: 

STEP  0:  Relax  (free)  all  integrality  restrictions  on  y^, 

prepare  for  storage  of  an  incumbent  solution,  set  the 
lower  bound  for  solution  value  to  +  •. 

REPEAT 

STEP  1:  Solve  ELUNLP  with  current  restrictions. 

STEP  2:  Improve  lower  bound  on  solution  value,  if  possible. 

STEP  3:  Round  solution  (heuristically)  to  satisfy  integer 
restrictions. 

STEP  4:  Compare  rounded  solution  with  the  incumbent,  replace 
incumbent  if  possible. 

STEP  5:  IF  current  solution  value  is  worse  than  current  lower 
bound,  go  to  STEP  7. 

STEP  6:  Fix  a  variable.  Heuristically  select  a  free  member  of 
yu  and  fix  it  at  a  bound.  Go  to  STEP  1. 

STEP  7:  Heuristically  select  STEP  8  or  STEP  9. 

STEP  8:  Reverse  a  variable.  Select  a  previously  fixed  variable 
and  reverse  it  to  its  opposite  bound,  go  to  STEP  1. 

If  no  fixed  variable  exists,  go  to  STEP  9. 


STEP  9:  Backtrack.  Select  a  reversed  variable  and  free  it.  Go 
to  STEP  1. 

UNTIL  TERMINATION. 

Note  that  the  heuristic  decision  rules  involve  the  method  of 
rounding  a  solution  and  the  selection  criteria  for  candidate  variables  to 
fix,  reverse,  and  free  (Figure  12).  Within  this  framework,  particular 
heuristics  yield  a  wide  variety  of  enumeration  strategies.  The  particular 
strategy  chosen  for  microcomputer  use: 

1.  rounds  the  current  restricted  solution  in  three  passes,  each  of 
which  selects  variables  from  a  class  defined  in  terms  of  e, 

where  0  £  e  £  .5  and  eu  £  yu  £  u  or  1  £  yu  £  el . 

Class  1:  nearly  integral  (0  £  e  £  .2) 

Class  2:  fractional  (.2  <  e£  .4) 

Class  3:  ambivalent  (.4  <  e£  .5) 

The  funding  heuristic  sequentially  exhausts  variables  from  each 
class  and  rounds  using  a  "minimal  regret  function,"  rounding  away 
from  the  worst  penalty. 


*  bound 


Fig.  12.  Primitive  Enumeration  Restrictions 
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2.  variables  are  selected  for  fixing  by  minimal  global  regret, 
and 

3.  variables  are  fixed,  reversed  and  freed  via  a  LIFO  List  operating 
exclusively  on  the  last  entry  in  the  list: 

Fix:  Push  fixed  variable  on  LIFO. 

Reverse:  Update  status  of  top  variable  on  LIFO. 

Free:  Pop  reversed  variable  from  LIFO,  mark  it  as 

freed  (Figure  12). 

Also  note  that  there  are  only  depth  and  value-motivated  fathoming 
rules  (decision  rules  for  reversal  or  backtracking);  feasibility  plays  no 
role  in  the  enumeration  except  via  the  value  of  the  Lagrangian  objective 
function. 

Finiteness  of  this  class  of  enumeration  methods  is  evident 
as  long  as  the  fix/reverse/backtrack  mechanism  cannot  yield  successive 
solutions  with  identical  restrictions. 

What  is  not  immediately  apparent  is  the  remarkable  effectiveness 
of  these  heuristics  with  the  elastic  enumeration  model.  Integer  solutions 
and  lower  bounds  of  excellent  quality  are  empirically  produced  quite 
early  in  the  enumeration  effort,  permitting  routine  early  termination  of 
the  search  based  on  an  optimality  tolerance  or  on  a  maximum  depth  (per¬ 
missible  number  of  fixed  variables  in  any  restriction);  tuning  of  the 
method  is  easily  accomplished  via  these  two  limits  and  the  elastic 
penalties  used  to  express  the  underlying  model. 
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or  Depth 


Sq  5  "relaxed"  solution,  V(S.)  solution  values 


Fig.  13.  Generic  Enumeration  Tree  (Fixed  Order) 


IV.  PACKAGE  DESCRIPTION 


The  mathematical  programming  package  described  here  provides  an 
integrated  repertoire  of  minimum-cost  network  flow  models  for  use  with  a 
microcomputer.  Various  interface  and  control  modules  are  also  described 
which  reduce  the  user's  workload  and  automate  the  solution  process. 
Extensive  facilities  for  data  file  management  and  provisions  for  incor¬ 
porating  custom  problem  generators  complete  the  package. 

A.  OVERVIEW 

1.  Package  Components 

Actually  a  suite  of  separate  programs  coordinated  by  a  master 
program,  the  package  is  automated  wherever  possible  and  completely 
interactive.  Designed  in  a  highly  structured  manner,  each  program  is 
modular,  relatively  standardized  and,  with  minor  modifications,  capable 
of  independent  operation.  Routines  common  to  more  than  one  program  or 
subject  to  frequent  modification  reside  in  the  system  library.  These 
features  simplify  modification  or  deletion  of  existing  modules,  addition 
of  new  programs,  or  even  the  incorporation  of  the  entire  package  (or  some 
subset)  into  a  larger  structure.  A  macro  view  of  the  package  is  given  in 
Figure  14. 

THE  APPLE  II  microcomputer  version  of  the  package,  exclusive  of 
data  files,  spans  two  disk  volumes  (removable  floppy  disks).  The  Pascal 
operating  system,  system  library,  master  program,  and  all  solution 
programs  reside  together  on  a  volume  that  is  always  on-line.  The  editor 
and  associated  preprocessor  programs  are  on  a  second  volume  that  can 
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I  I 

|  SYSTEM*  | 
I  LIBRARY  | 


I  I 

|  MASTER  | 
|  PROGRAM  | 


PREPROCESSORS 

*  ACCESSIBLE  TO  ALL  PROGRAMS 


Fig.  14.  Package  Block  Diagram 
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be  taken  off-line  after  the  desired  program  is  loaded  into  main  memory. 
Oata  files  occupy  additional  volumes  present  only  for  read/write  opera¬ 
tions.  This  partition  is  dictated  by  the  size  of  the  package  code  and 
has  been  chosen  to  minimize  disk  manipulation  requirements  in  a  two-drive 
system  (the  minimum  practical  configuration).  Package  operation  is  thus 
able  to  proceed  with  volume  exchange  required  only  for  data  file  access 
and  editing  procedures. 

2.  Data  Files 
a.  Structure 

Although  several  quite  distinct  models  are  supported,  a 
common  format  is  possible  through  the  use  of  Pascal's  facility  of  variant 
records.  A  Pascal  record  is  a  compound  structure  of  arbitrary  types  of 
data  which,  when  composed  of  types  with  partly  identical  components,  is 
termed  a  variant  record.  A  portion  of  the  record  is  the  same  for  all 
occurences  while  the  remainder  (the  variant  part)  may  differ  depending  on 
the  value  of  an  indicator  variable  (also  part  of  the  record).  Wirth 
[Ref.  34]  gives  an  excellent  description  of  such  data  structures  and 
their  employment.  This  allows  use  of  a  single  data  structure  for  inter¬ 
program  transfer  of  information  with  specific  portions  of  the  package 
extracting  only  those  items  they  actually  need.  Each  program  then 
converts  this  input  data  into  the  required  internal  data  representation. 
Solution  technique  selection  can  thus  be  accomplished  without  user 
intervention  as  the  data  record  contains  enough  information  to  determine 
problem  class  and  a  record  can  be  accessed  by  any  portion  of  the  package. 
The  use  of  the  same  type  record  for  all  problems  provides  streamlined 
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data  access  procedures  without  sacrificing  the  efficiency  of  custom 
internal  data  structures  for  each  program. 

A  typical  data  file  is  a  collection  of  randomly  accessable 
records,  each  of  which  has  three  possible  structures:  (1)  header,  (2) 
arc,  and  (3)  node.  The  Pascal  interpretation  of  this  scheme  is  described 
in  Appendix  A.  Although  only  one  header  record  is  allowed  per  file,  any 
number  of  arc  and  node  records  (subject  to  free  disk  volume  space)  can  be 
contained  in  a  single  file  without  regard  to  order.  All  arcs  of  the 
model  must  be  explicitly  represented  in  the  data  file;  however,  only 
those  nodes  with  attributes  not  equal  to  default  values  must  be  included. 
Solution  programs  assign  appropriate  default  values  to  all  nodes  and  then 
process  the  input  file  making  note  of  the  nodes  that  deviate  from  default 
settings.  A  disk  volume  dedicated  to  one  file  can  accormrodate  approxi¬ 
mately  three  thousand  such  records  (APPLE  II). 

Constructed  by  the  EDITOR  during  file  creation,  the  header 
record  describes  the  problem  represented  by  the  file.  Problem  name, 
problem  type,  number  of  nodes,  number  of  arcs,  and  access  history  are 
included  in  this  record. 

Each  arc  record  contains  the  arc  name,  source  and  destination 
nodes,  bounds  on  flow,  initial  flow,  and  cost  information.  Depending  on 
the  type  of  arc  represented,  the  cost  data  varies.  For  linear  cost 
functions,  only  the  cost  per  unit  flow  and  the  "1-u"  status  are  needed. 
Nonlinear  cost  functions  require  specification  of  the  function  identifica¬ 
tion  (assigned  number  in  the  function  library)  and  coefficient  structure. 

Node  records  specify  name,  identification  number  in  the 
model,  type,  flow  requirements,  flow  range,  and  penalty  for  range  violation. 
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The  identification  number  must  be  unique,  for  a  given  model  (enforced  by 
the  EDITOR  on  file  creation  and  update).  Flow  requirements  determine 
the  node  type  as  follows:  zero  f low--transshipment,  positive  flow- 
supply,  and  negative  flow-demand. 

Internal  data  representations  will  be  presented  as  the 
various  programs  are  discussed. 

b.  File  Naming  Conventions 

In  order  to  reduce  operator  workload,  a  menu-driven  data 
file  selection  technique  is  employed.  Each  data  file  name,  regardless 
of  the  type  model  it  represents,  contains  the  suffix  ".net."  Whenever  a 
file  is  to  be  accessed,  a  directory  for  the  on-line  volumes  is  displayed, 
files  with  this  suffix  are  identified,  and  a  single  keystroke  selects  the 
desired  file.  This  approach  ensures  that  only  files  compatible  with  the 
package  are  accessed  and  greatly  simplifies  their  retrieval  (from  the 
user's  point  of  view).  The  EDITOR  program  appends  the  required  suffix  to 
user-designated  file  names  upon  file  creation. 

3.  User  Friendly  Features 

Menu-driven  displays  and  single  keystroke  option  selection 
make  the  package  easy  to  use.  On  input,  all  user  entries  are  examined 
for  errors  in  the  sense  of  range,  type  matching  (numeric  vs.  alpha),  and 
context  as  applicable. 

Two  general  menu  formats  are  standard  throughout  the  package. 

The  first  is  used  for  selecting  courses  of  action.  In  such  a  situation, 
one  of  the  displayed  items  must  be  selected  before  the  program  continues. 
The  option  selection  menu  from  the  EDITOR  program  serves  as  an  example  of 
this  process  and  is  shown  in  Figure  15.  A  common  decision  point  in 


EDITOR  PROGRAM 


OPTIONS 

####### 

A(Uer  an  existing  file 
B(rowse  through  file 
C(reate  a  new  file 
R(emove  a  file  [ permanently] 

Transfer  a  file 

<ESCAPE>  EDITOR  program  and  return  to  MASTER  program 
OPTION  DESIRED 

— c  ]— 


Fig.  15.  Typical  Option  Selection  Menu 

optimization  programs  concerns  the  setting  of  program  parameters.  The 
second  menu  format  deals  with  this  by  displaying  the  default  values  or 
choices  as  appropriate.  To  change  a  particular  value,  the  user  enters 
the  menu-designated  symbol  associated  with  the  parameter  in  question.  If 
the  parameter  is  an  ON/OFF  choice,  the  change  is  made  when  the  symbol  is 
entered;  for  numerical  quantities,  the  user  is  prompted  for  the  new 
value.  Menu  updating  occurs  automatically  until  the  user  indicates  that 
all  values  are  correct.  Figure  16  illustrates  a  typical  parameter 
control  menu.  Choices  from  either  format  may  invoke  further  sub-menus. 
Only  selections  contained  in  the  current  menu  are  admitted:  when  the 
user  enters  a  choice  symbol  not  depicted  by  the  menu  in  use,  an  appro¬ 
priate  error  message  is  issued. 


ELASTICNET 

OPTIONS 

####### 

C(omplete  printout  --> 

ON 

D(etailed  printout  --> 

OFF 

F(ile  output  — > 

OFF 

H(ard  copy  --> 

OFF 

Internal  array  dump  — > 

OFF 

M(odel  alteration  — > 

OFF 

Are  options  satisfactory  ? 

Y(es  or  N(o  —  > 

Fig.  16.  Typical  Parameter  Control  Menu 


B.  CONTROL  AND  UTILITY  PROGRAM  DESCRIPTIONS 
1.  Master  Program 

The  master  program  orchestrates  package  operation  by  passing 
control  to  appropriate  programs  as  requested  by  the  user  (either  directly 
or  indirectly).  This  is  accomplished  with  two  levels  of  command  internal 
to  the  master  program  as  shown  in  Figure  17.  The  outer  command  level 
allows  the  choices  of  fl)  data  file  manipulation,  (2)  problem  solution, 
and  (3)  package  exit.  The  first  and  third  choices  transfer  control  to 
the  EDITOR  program  and  the  Pascal  operating  system  respectively.  A 
problem  solution  request  activates  the  solution  command  level  and  package 
flow  proceeds  as  depicted  in  Figure  18. 

Problem  solution  begins  with  user  selection  of  the  input  file. 


Once  the  data  file  has  been  identified,  the  course  of  action  is  determined 
by  problem  size  (number  of  nodes,  number  of  arcs)  and  type  as  indicated 


Fig.  17.  Outer  Command  Logic 
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*  indicates  user  interaction  required  (otherwise  automatic) 
f  indicates  elastic  &  "o-u" 

£  indicates  warning  message  issued 

Fig. 18.  Solution  Process  logic 
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by  the  file  header  record.  Linear  problem  files  are  those  with  all 
linear  arc  cost  functions  to  include  problems  with  explicit  elastic 
ranged  or  "1-u"  constraints.  Nonlinear  elastic  or  "l-u“  problems  are  not 
supported  and  are  treated  as  ordinary  nonlinear  models. 

All  linear  problems,  size  permitting,  are  routed  to  ELASTICNET 
(the  combination  ENET/UNET  program),  otherwise  GNET  is  invoked  (again 
size  permitting).  Nonlinear  problems  are  first  examined  to  ensure  that 
their  size  is  commensurate  with  NLPNET  capabilities  and  that  all  required 
cost  functions  are  resident  in  the  system  library.  If  either  one  of 
these  tests  fails,  the  user  is  given  the  option  to  linearize  and  use 
only  GNET  for  an  approximate  solution  (assuming  size  is  within  GNET 
limits).  Once  it  has  been  decided  that  NLPNET  can  safely  be  called, 
initial  flow  feasibility  is  determined  by  a  straightforward  node  flow 
conservation  calculation.  Feasible  problems  go  directly  to  NLPNET  for 
solution,  however,  infeasible  problems  are  temporarily  linearized  and  a 
feasible  starting  point  obtained  by  GNET  prior  to  being  sent  to  NLPNET. 

In  the  event  that  GNET  determines  that  no  feasible  solution  exists,  the 
solution  process  is  terminated.  Error  messages  are  issued  whenever  the 
outcome  of  the  process  differs  from  that  expected  for  the  type  file  in 
question. 

A  variety  of  elementary  linearization  options  are  offered  as 
noted  in  Figure  19.  The  arc  cost  functions  are  evaluated  at  user  desig¬ 
nated  points  and  combined  as  specified  to  obtain  point  or  interval 
approximations  which  then  serve  as  cost  coefficients  for  GNET.  These 
costs  are  discarded  after  use  by  GNET.  It  should  be  noted  that  this  is 
not  a  piecewise  linearization  and  is  not  intended  to  be  a  complete 
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1.  Point  gradient  eval  @  LCC*  of  bounds  (t*l  +  (l-t)*u). 

2.  Point  gradient  eval  @  midpoint  of  bounds  (t  =  0.5). 

3.  Point  gradient  eval  @  LCC  of  lowerbound  and  present  flow. 

4.  Point  gradient  eval  @  LCC  of  present  flow  and  upperbound. 

5.  Interval  linearization  from  lowerbound  to  specific  point. 

6.  Interval  linearization  over  entire  bound  interval. 

7.  Assign  zero  cost  to  all  arcs. 

*  LCC  =  Linear  Convex  Combination 

OPTION  DESIREO 
— [ 


Fig.  19.  Linearization  Options 

solution  technique.  GNET  can  support  piecewise  linearizations  if  the 
model  is  explicitly  described  by  addition  of  the  appropriate  arcs.  The 
user  can  request  that  all  nonlinear  problems  be  linearized,  regardless  of 
feasibility  status,  by  overriding  a  default  parameter. 

The  program  calls  depicted  in  Figure  18  are  automatically  gener¬ 
ated  by  the  solution  control  module  of  the  MASTER  program.  Information 
such  as  data  file  names  and  package  control  instructions  are  passed  to 
activated  programs  through  the  Pascal  operating  system.  Once  control  is 
passed  to  a  program,  this  information  is  retrieved  and  the  subject  file 
is  verified  as  still  being  on-line.  In  the  event  that  the  volume  contain¬ 
ing  the  data  file  has  been  moved  off-line,  the  program  issues  an  error 
message  and  awaits  user  action. 

2.  EDITOR  Program 

The  EDITOR  program  allows  the  user  to  create,  alter,  transfer,  and 
examine  data  files  as  indicated  in  Figure  20.  Since  this  portion  of  the 
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Fig.  20.  Editor  Block  Diagram 


package  requires  the  most  user  participation,  every  attempt  has  been  made 
to  make  the  program  as  user-friendly  as  possible.  In  addition  to  employ¬ 
ing  the  standard  interactive  techniques,  the  EDITOR  program  input  format 
can  be  reconfigured  to  conform  to  the  user's  needs.  This  redefinition  of 
the  input  field  order  can  be  effected  at  any  time  during  the  editing 
process.  For  example,  one  file  can  be  created  utilizing  the  default 
input  order,  reconfiguration  performed,  and  a  second  file  altered  with  a 
completely  different  input  field  order.  Instead  of  reordering  the  raw 
data,  one  merely  changes  the  order  in  which  the  program  expects  data 
fields  to  be  received.  Additionally,  fields  can  be  omitted  entirely  from 
the  input  format  and  default  or  user-definable  constant  values  assigned. 
Input  field  order  for  arcs  and  nodes  may  differ. 

a.  File  Creation 

Files  can  be  created  from  text  (human  readable)  files, 
keyboard  input,  or  through  the  use  of  preprocessor  programs.  For  all 
modes  of  file  creation,  the  program  keeps  track  of  the  filetype  by  typing 
the  arcs  as  they  are  received  (i.e.,  one  nonlinear  arc  changes  a  linear 
file  to  nonlinear).  Input  from  text  files  and  the  keyboard  is  examined 
for  the  following  inconsistencies: 

1.  Prohibited  node  numbers  (nonpositive  or  greater  than  number 
of  nodes. 

2.  Inconsistent  upper  and  lower  limits  on  quantities  (lower 
limit  not  less  than  or  equal  to  upper  limit). 

3.  Assignment  of  quantity  values  outside  of  established 
limits. 
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b.  File  Alteration 

Records  can  be  added,  deleted  and  their  individual  data 
fields  updated.  Addition  of  records  is  allowed  for  both  node  and  arc 
records;  however,  only  arcs  may  be  deleted.  This  restriction  is  imposed 
because  no  connectivity  analysis  is  performed  by  the  EDITOR  on  the 
network  representation  resulting  from  alterations.  Deletion  of  arcs  can 
isolate  portions  of  the  network,  but  the  solution  programs  acconmodate 
this  situation  effectively.  Removing  nodes  requires  an  exhaustive  search 
of  the  data  file  to  detect  arcs  that  must  also  be  deleted.  This  capabil¬ 
ity  is  not  presently  supported  by  the  EDITOR. 

c.  File  Transfer 

Data  files  can  be  transferred  to  the  console,  printer, 
another  data  file,  a  text  file,  or  a  remote  computer.  The  data  field 
format  in  effect  at  transfer  governs  the  order  of  field  transmission  so 
the  package  can  be  linked  to  data  bases  with  diverse  format  requirements. 
File  transfers  to  and  from  other  computers  are  limited  to  text  files 
only,  thus  the  ability  to  convert  between  data  files  and  text  files  has 
been  provided.  A  stand-alone  file  transfer  program  [Ref.  35]  has  been 
used  for  such  communication  with  great  success. 

3.  Preprocessor  Programs 

These  programs  are  created  for  specific  models  to  assist  in 
data  base  reduction.  For  example,  if  data  on  a  certain  network  exists  in 
raw  form  such  as  physical  measurements  on  the  actual  entities  represented 
by  the  nodes  and  arcs  of  the  model,  the  user  can  write  a  custom  program 
to  translate  such  information  to  a  form  suitable  for  package  use  and  then 
execute  the  conversion  directly  from  the  EDITOR  program.  A  catalog  of 
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existing  preprocessor  programs  is  maintained  by  the  EDITOR  program  and 
presented  to  the  user  upon  request.  After  selecting  the  new  data  file 
name,  control  is  passed  to  the  desired  preprocessor  and  returned  to  the 
EDITOR  upon  completion  of  preprocessing. 

Facilities  are  provided  by  the  EDITOR  to  maintain  the  preprocessor 
catalog  (a  data  file)  and  manually  activate  uncataloged  preprocessors 
(which  are  then  automatically  added  to  the  catalog).  The  preprocessors 
are  stand-alone  programs  not  constrained  to  reside  on  one  of  the  package 
volumes  with  no  (package  imposed)  limits  on  their  number  or  individual 
size. 

C.  SOLUTION  PROGRAMS 

Each  program  is  partitioned  into  four  segments  as  shown  in  Figure  21. 
The  main  segment  (of  the  particular  solution  program)  calls  subordinate 
segments  into  memory  sequentially,  thereby  maximizing  memory  available  for 
data.  The  input,  initialization,  and  report  segments  contain  all  proce¬ 
dures  that  communicate  with  either  the  user  or  other  portions  of  the 
package.  Solution  segments  are  therefore  minimal  representations  of 
their  respective  algorithms  and  data  structures,  and  are  devoid  of  any 
nonessential  information. 


1 

1 

INPUT* 

MAIN 

1 

INITIALIZATION* 

SEGMENT 

1 

1 

SOLUTION 

i 

! 

REPORT* 

*User-friendly  structures  used 


Fig.  21.  Generic  Solution  Program  Segment  Map 
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The  report  segment  automatically  updates  the  data  file  with  the 
new  solution  and  displays  the  results  of  the  optimization  on  user- 
designated  peripherals.  Detailed  reports  of  the  solution  process  may  be 
provided  upon  user  request.  Subsequent  program  calls  are  controlled  by 
logic  contained  in  the  report  segment. 

1.  GNET 

This  program  is  constructed  so  that  both  linear  and  nonlinear 
problems  are  processed  in  essentially  the  same  manner.  Two  minor  dif¬ 
ferences  are  accommodated  by  a  flag  that  keys  on  the  problem-type  field 
of  the  header  record.  First,  nonlinear  files  sent  to  GNET  for  solution 
have  the  linearized  arc  costs  placed  in  the  initial-flow  field  of  the  arc 
data  record  as  GNET  does  not  use  this  field  for  linear  problem  input 
data.  Additionally,  on  obtaining  a  feasible  solution,  control  is  passed 
to  NLPNET  vice  the  MASTER  program  when  GNET  terminates  from  a  nonlinear 
problem. 

a.  Input  Requirements 

Upon  activation,  GNET  extracts  the  information  listed  in 
Figure  22  from  the  input  data  file.  Certain  compacting  operations  take 
place,  as  each  arc  record  is  processed,  to  reduce  storage  requirements. 
The  arc  flow  range  is  stored  as  one  number,  the  upper  bound,  after 
translation  by  the  lower  bound.  Also,  the  destination  node  list  is 
maintained  with  all  arcs  having  the  same  destination  node  stored  in 
contiguous  locations.  This  allows  the  list  to  be  a  node-length  array 
instead  of  an  arc-length  array.  The  arc  costs  and  source  nodes  are 
read  directly  into  arc-length  arrays. 


ARCS: 


Source  node  number 
Destination  node  number 
Lower bound  on  flow 
Upperbound  on  flow 
Cost  information 

NODES: 


Net  flow 


Fig.  22.  GNET  Input  Data  Requirements 

Node  information,  other  than  net  flow,  is  ignored  as  GNET 
does  not  support  ranged  constraints.  In  the  event  that  total  supply  does 
not  equal  total  demand,  the  program  will  terminate  abnormally.  GNET 
operates  with  integer  arithmetic;  therefore  all  input  values  not  explic¬ 
itly  integer  are  truncated  upon  receipt, 
b.  Internal  Data  Structure 

GNET  introduces  an  artificial  node  called  the  "root"  to 
complete  the  basis  tree.  Three  node  length  arrays  are  used  to  represent 
the  basis  tree  and  provide  a  mechanism  to  easily  traverse  the  tree.  The 
predecessor  array  defines  the  next  node  on  the  path  from  a  given  node  to 
the  root.  For  example,  predecessor  (i)  contains  the  node  number  of  the 
predecessor  of  the  ith  node.  The  depth  array  stores,  for  each  node, 
the  number  of  nodes  on  the  path  to  the  root.  Finally,  the  traversal 
array  provides  for  each  node  the  next  node  to  evaluate  during  forward- 
substitution  of  the  basis.  This  array  gives  a  sequence  of  nodes  which, 
combined  with  the  respective  node  predecessors,  define  an  upper  tri angu¬ 
lation  of  the  basic  arcs.  The  ith  element  of  this  array  is  the  next 
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node  that  would  be  visited  from  node  i  in  a  classical  preorder  tra¬ 
versal  .  Arrays  to  store  flow  on  basic  arcs  and  the  simplex  dual  variables 
complete  the  data  structure.  Figure  23  illustrates  the  basis  representa¬ 
tion  and  Figure  24  describes  the  various  arrays. 

These  arrays  are  used  by  GNET  to  perform  the  simplex  opera¬ 
tions  with  elementary  algebraic  and  logic  operations.  A  detailed 
description  of  this  data  structure  and  its  employment  is  given  by 
Bradley,  Brown,  and  Graves  [Ref.  3]. 
c.  Solution  Segment 

An  all-artificial  starting  technique  is  used  to  begin  the 
solution  process.  The  initial  basis  consists  of  artificial  arcs  between 
each  node  and  the  root  node.  The  flow  on  these  arcs  is  set  equal  to  the 
demand  or  supply  requirements  of  the  respective  node.  Arcs  are  oriented 
from  the  root  to  demand  nodes  and  to  the  root  from  supply  nodes.  Assign¬ 
ing  a  relatively  large  cost  to  each  artificial  arc  ensures  that  a  feasible 
basis  will  contain  no  artificial  arcs  at  positive  flow.  Artificial  arcs 
present  at  the  conclusion  of  the  optimization  with  non-zero  flow  are 
explicitly  identified. 

2.  NLPNET 

a.  Input  Oata  Requirements 

NLPNET  requires  the  input  data  illustrated  in  Figure  25. 

The  source  node,  destination  node,  arc  flow  bounds,  and  initial  flow  are 
read  directly  into  arc-length  arrays.  Since,  in  a  nonlinear  problem,  all 
arcs  potentially  can  have  non-zero  flow  at  optimality,  the  flow  array 

"*As  defined  in  the  computer  science  literature  [e . g . ,  Ref.  36]. 
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Node: 


2  3  4  5  6  7 


Predecessor: 
Traversal  : 
Depth  : 


I 

3  5  -5  3  7  -5  A  | 

4  6  1  2  3  7  5  1 

3  2  2  3  1  2  0! 


— ->  Traversal 


Fig.  23.  GNET  Basis  Representation 
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Problem  Definition  (All  Arc-Length  Arrays  Except  as  Noted) 


H 


T 

C 

CP 


(  )  -  A  node-length  array  of  entry  points  into  arc-length  arrays 
for  each  head  node  (e.g.,  all  arcs  directed  toward  head 
node  i  are  found  in  locations  H(i ),...,  H(i  +  1)  -  1 
of  arc  length  arrays). 

(  )  -  Tail  node  indices  ( i . e . ,  nodes  which  arcs  are 
directed  away  from). 

(  )  -  Cost  per  unit  flow. 

(  )  -  Upper  bound  (capacity)  on  flow.  (The  sign  bit  is  used 
to  indicate  a  reflected  arc.) 


Basis  Representation  (All  Node-length  Arrays) 


P  (  ) 


IT  (  ) 


D  (  ) 


The  predecessor  of  each  node  on  its  backpath  to  the  root 
node.  The  orientation  of  the  basic  arc  connecting 
node  i  to  its  predecessor  is  indicated  by  the  sign 
of  P(i),  a  negative  value  indicating  an  arc 
(i,  -  P( i ) ) »  and  a  positive  value  signifying  an 
arc  (P(i),  i). 

Preorder  traversal  successors.  IT(i)  gives  the 
node  number  of  the  next  node  to  visit  in  preorder 
after  node  i.  With  P(  ),  IT(  )  defines  a  basis 
triangulation  and  can  be  used  for  substitution 
solution. 

Depths  of  the  nodes  in  the  current  basis.  D(i) 
gives  the  length  of  the  backpath  from  node  i  to 
the  root  node  (used  in  column  generation). 


Solution  Representation  (All  Node- length  Arrays) 

X  (  ),-  Current  flow  of  each  basic  arc,  and  capacity  minus 
CPX  (  )  current  flow  of  each  basic  arc  (i.e.,  if 

P(i )  3  J  <  0,  flow  is  X(i )) . 


Fig.  24.  GNET  Data  Structure 


AD- A 109  699  NAVAL  POSTGRADUATE  SCHOOL  MONTEREY  CA 

A  MICROCOMPUTER-BASED  NETWORK  OPTIMIZATION  PACKAGE. (U) 
SEP  81  R  H  DUFF 


1.0  s 


1.1 


.25 


IM  IM 

“  Itt 

*-  HIM 

■  18 

1.6 


1.4 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BURTAU  Of  STANDARDS  196T  A 


ARCS: 


Source  node  number 
Destination  node  number 
Lowerbound  on  flow 
Initial  flow 
Upperbound  on  flow 
Cost  function  type 
Cost  function  coefficients 

NODES: 


Net  flow 


Fig.  25.  NLPNET  Input  Data  Requirements 

must  be  an  arc-length  structure.  The  bounds  and  the  node  pairs 
associated  with  each  arc  are  explicitly  maintained,  in  contrast  to 
the  GNET  scheme  where  non-basic  arcs  are  coded  with  a  single  bit  to 
indicate  status  (at  upper  or  lower  bounds). 

Nonlinear  cost  function  definitions  reside  in  the  SYSTEM 
LIBRARY  where  they  can  easily  be  modified  without  requiring  changes  to 
any  of  the  programs  comprising  the  package.  Each  function  is  assigned 
a  unique  number  upon  inclusion  in  the  library,  so  determination  of  a 
function's  presence  is  a  trivial  matter.  This  is  accomplished  by  checking 
function  identification  numbers  against  a  master  list  located  in  the 
SYSTEM  LIBRARY.  Also  required  in  the  function  definition  are  analytical 
expressions  for  the  gradient  and  hessian.  Functions  can  be  defined  using 
up  to  three  coefficient  terms;  if  more  than  three  such  coefficients  are 
needed,  then  a  small  coding  change  in  NLPNET  is  required. 

Input  node  information  is  the  same  as  for  GNET  (net  flow 
only  for  non-default  nodes).  As  in  GNET,  total  supply  and  demand 
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conservation  must  be  achieved  or  the  program  terminates  at  the  input 
stage. 

b.  Internal  Data  Structure 

NLPNET  also  uses  a  root  node  to  complete  the  basis.  The 
basis  tree  is  depicted  by  the  familiar  predecessor,  depth,  and  traversal 
arrays  which  are  functionally  equivalent  to  the  corresponding  6NET 
structures.  A  fourth  array,  a  reverse-traversal  array,  allows  mobility 
opposite  to  that  provided  by  the  traversal  array.  This  array  is  the 
inverse  of  the  preorder  traversal  and  is  used  for  back-substitution 
during  the  calculation  of  the  basic  variable  search  direction  induced  by 
a  superbasic  direction.  The  basis  representation  is  described  by 
Figure  26.  Arrays  to  maintain  the  variable  partition,  gradient  vector, 
hessian  diagonal  vector,  search  direction,  and  dual  variable  estimates 
complete  the  data  structure.  These  arrays  allow  NLPNET  to  perform  primal 
simplex  procedures  directly  on  the  basis  representation  in  the  spirit  of 
GNET.  A  summary  of  the  data  structure  is  shown  in  Figure  27.  A  complete 
description  of  the  arrays  and  t.ieir  use  is  given  by  Dembo  [Ref.  19]. 

c.  User  Definable  Parameters 

Real  arithmetic  is  used  extensively  throughout  NLPNET  so 
various  tolerances  are  necessary.  Additionally,  the  solution  process 
employs  several  parameters  that  control  the  performance,  convergence,  and 
accuracy  characteristics  of  the  solution  process.  Default  settings  for 
these  values  are  built  into  the  program,  however,  the  user  may  redefine 
these  settings  (within  established  limits)  if  desired.  Figure  28  is  the 
menu  that  the  user  receives  when  the  default  settings  are  to  be  modified. 
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Node  i 

1 

2 

3 

4 

5 

6 

Predecessor*  S 

3 

4 

5 

-5 

6 

6 

Traversal  1 

6 

3 

1 

2 

4 

5 

Inverse  Traversal  | 

3 

4 

2 

5 

6 

1 

Depth  | 

1. 

3 

3 

2 

2 

1 

0 

— >  Traversal  Path 
••••>  Inverse  Traversal  Path 

♦Negative  I  =  predecessor  (J)  implies  that  the  arc  is  orientated 
from  node  I  to  node  J. 


Fig.  26.  NLPNET  Basis  Representation 
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Problem  Definition  Arrays 


IFROM  (  )  - 
ITO  (  )  - 
LOWER  (  )  - 
UPPER  (  )  - 
ITYPE  (  )  - 

COEF  (  )  - 
IPTC  (  )  - 

RHS  (  )  - 


Source  nodes  of  arcs. 

Destination  nodes  of  arcs. 

Lower  bounds  on  arc  flow. 

Upper  bounds  on  arc  flow. 

Library  identification  numbers  of  arc  objective 
functions. 

Coefficients  of  objective  function  (0-3  per  arc). 
Pointers  to  index  of  coefficient  array  contains  first 
coefficient  of  each  objective  function. 

Node  net  flow  requirements 


Variable  Partition 


LB  (  )  - 
LS  (  )  - 
LN  (  )  - 
ISBEST (  )  - 


Basic  arc  indices. 

Superbasic  arc  indices. 

Nonbasic  arc  indices. 

Best  superbasic  variable  to  replace  a  given  basic 
variable. 


Basis  Representation 

LTHRO  (  )  -  Traversal  array  (thread). 

LPRED  (  )  -  Predecessor  array. 

LRTHRD(  )  -  Inverse  of  traversal  array  (reverse  thread). 
NDEPTH(  )  -  Depth  array. 


Solution  Representation 


FLOW 

F 

G 

H 

P 

W 


(  )  -  Arc  flow. 

(  )  -  Arc  contribution  to  objective  function  (at  FLOW). 
(  )  -  Gradient  vector. 

(  )  -  Hessian  diagonal  vector.  (Hessian  is  a  diagonal 
matrix. ) 

(  )  -  Search  direction  vector. 

(  )  -  Dual  variable  estimates. 


Fig.  27.  NLPNET  Data  Structure 


CURRENT  PARAMETER  VALUES 


GENERAL  TOLERANCES 


1.  Alpha  =  O.OOOl 

2.  Flow  =  0.0001 

3.  Gradient  =  0.001 

4.  Price-out  =  0.0001 

LINESEARCH  TOLERANCES 

5.  Eps  *  0.001 

6.  T  =  0.02 

7.  Eta  =  0.1 

8.  Epsd  =  0.001 

OTHER 

9.  Forcing  ftn  =1.0 

10.  Max  Iterations  =  50 


Parameters  OK?  :  Y(es  or  symbol  to  change  — > 


Fig.  28.  NLPNET  Oefault  Modification  Menu 
d.  Reports 

The  program  prints  iteration  reports  and  final  solution 
reports  without  user  intervention.  Either  of  these  reports  may  be 
selectively  disabled.  In  any  event,  the  input  data  file's  initial  flow 
fields  will  always  be  updated  upon  normal  program  termination.  Appendi 
B  describes  these  reports. 


3.  ELASTICNET 

This  program  combines  both  ENET  and  UNET  into  one  composite 
solution  module  that  is  resident  in  memory  for  the  duration  of  the 
optimization  process.  In  this  manner,  costly  disk  access  operations 
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resulting  from  UNET  communication  with  ENET  during  "l-u"  problem  solution 
are  eliminated.  Of  course,  such  a  scheme  requires  more  memory  resources 
for  code  storage  than  would  be  the  case  with  a  partitioned  format;  however, 
UNET  is  a  very  concise  code  and  the  storage  reduction  is  minimal. 

The  selection  of  the  appropriate  solution  mechanism  is  auto¬ 
matically  controlled  by  logic  internal  to  ELASTICNET  activated  by  the 
presence  of  "l-u"  arcs  in  the  input  data  file.  Invocation  of  UNET 
in  response  to  such  arcs  can  be  suppressed  from  the  program  option  menu 
and  only  ENET  utilized  to  process  a  file  with  "l-u"  arcs  (giving  only  a 
relaxed  solution).  This  option  allows  a  single  file  to  represent  both  a 
free  and  "l-u"  model. 

a.  Input  Data  Requirements 

ELASTICNET  exercises  virtually  the  entire  input  data  struc¬ 
ture  and  extracts  the  data  given  in  Figure  29  from  the  data  file.  For 
arcs,  the  source  node,  destination  node,  and  unit  cost  are  read  into 
arrays  identical  to  the  GNET  counterparts  already  described.  ELASTICNET, 
unlike  GNET,  maintains  explicit  representations  for  each  arc's  upper 
and  lower  bounds  on  flow.  This  simplifies  the  program  and  results  in  a 
miniscule  increase  in  storage  requirements.  Finally,  the  "l-u"  status  of 
each  arc  is  maintained  in  an  arc-length  array  with  zeros  indicating  free 
arcs  and  ones  marking  "l-u"  arcs. 

Input  node  information  consists  of  five  numbers  for  each 
node:  net  flow  at  the  node,  upper  and  lower  range  on  flow,  and  penalties 
for  upper  and  lower  range  violation.  In  the  absence  of  input  information 
for  ranges  and  penalties,  default  values  are  assigned  by  the  input 
module.  As  with  the  other  programs,  all  supply  and  demand  nodes  must 
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ARCS: 


Source  node  number 
Destination  node  number 
Lower bound  on  flow 
Upperbound  on  flow 
Cost  information 
"l-u"  status 


NODES: 


Net  flow 

Lower  range  on  flow 
Upper  range  on  flow 
Lower  penalty  for  range  violation 
Upper  penalty  for  range  violation 


Fig.  29.  ELASTICNET  Input  Data  Requirements 

explicitly  be  present  in  the  input  file.  Additionally,  any  transshipment 
nodes  for  which  ranges  or  penalties  differ  from  default  values  must  be 
included.  The  magnitude  of  the  default  settings  is  user  definable  from 
the  option  menu.  Flow  conservation  (supply  »  demand)  is  not  required 
for  the  elastic  model  accommodated  by  ELASTICNET  and  therefore  only  a 
warning  message  is  issued  when  supply  does  not  equal  demand, 
b.  Internal  Data  Structure 

ELASTICNET  maintains  a  basis  representation  (Figure  30) 
very  similar  to  the  GNET  structure.  Additional  arrays  are  required  to 
control  the  UNET  enumeration  process  and  record  incumbent  solutions. 

These  structures  are  described  in  Figures  31  and  32. 
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Node :  i 

1 

2 

3 

4 

5 

6 

7 

Predecessor : 

1 

3 

5 

-5 

3 

-7 

-5 

a 

Traversal  : 

! 

1 

2 

7 

4 

3 

6 

5 

Depth  : 

1 

0 

0 

3 

0 

2 

0 

1 

Aggregate  : 

1 

1 

C13 

C25 

2 

C34 

2 

C56 

0 

Fig.  30.  ENET  Basis  Representation 
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4  _ _ 


T 


rvr 


Problem  Definition 


Arcs: 

(All  Arc-Length  Arrays  Except  as  Noted) 

H 

( 

)  - 

A  node-length  array  of  entry  points  by  head 

node  ii 

arc  lists  (sign  bit  indicates  S.  reflected), 
Tail  node  indices  (sign  bit  indicates  fixed 

T 

( 

)  - 

arc) . 

C 

( 

)  - 

Cost  per  unit  flow. 

BL 

( 

)  - 

Lower  bounds  on  flow. 

BU 

( 

)  - 

Upper  bounds  on  flow  (sign  bit  indicates  reflected 

Nodes 

: 

(All  Node-Length  Arrays) 

RL 

( 

)  - 

Lower  ranges. 

RU 

( 

)  - 

Upper  ranges. 

DL 

( 

)  - 

Penalties  for  lower  range  violation. 

DU 

( 

)  - 

Penalties  for  upper  range  violation. 

Basis  Representation  (Similar  to  GNET) 

P  (  )  -  Predecessor  array  (sign  bit  indicates  basic  arc 
orientation) . 

IT  (  )  -  Traversal  array. 

D  (  )  -  Depth  array. 

A  (  )  -  Aggregated  successor  array  (for  an  aggregated  node, 

the  cost  is  stored  for  its  basic  arc  predecessor; 
for  a  disaggregated  node,  the  number  of  aggregated 
successor  nodes/basic  arcs  is  stored)  [Ref.  3,  p.  28]. 

Solution  Representation  (All  Node-Length  Arrays) 


X 

( 

),- 

Arrays  with  current  flow  of  each  basic  arc,  and 

8UX 

( 

) 

capacity  minus  current  flow  of  each  basic  arc 
(i.e.,  P(i)  *  j  <  0,  flow  is  X( i ) ) . 

U 

( 

)  - 

Dual  variables. 

Fig.  31.  ENET  Data  Structure 
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Enumeration  Control 

IFIX(  )  -  A  last-in-first-out  (LIFO)  structure  that  records  the 
arcs  currently  fixed  at  a  bound.  IFIX(i)  is  the  arc 
number  in  the  current  enumeration  of  the  itn 
restriction. 

LGB  {  )  -  Maintains  (3-bit)  "1-u"  status  of  each  arc.  On  input. 

"1-u"  arcs  are  assigned  a  LGB  value  of  one  while  all 
other  arcs  are  given  a  value  of  zero.  During  enumer¬ 
ation,  nonzero  values  of  LGB  indicate  bounds  at  which 
arcs  are  fixed  or  reversed.  Fixed  arcs  have  negative 
T  (  )  while  free  arcs  have  positive  T  (  )  values. 

Incumbent  Solution  Record 

IPS  (  )  -  The  best  primal  solution  encountered. 

IDS  (  )  -  The  dual  solution  corresponding  to  the  best  primal 
solution. 

IB  (  )  -  Initial  right  hand  side  of  constraints. 


Fig.  32.  UNET  Data  Structure 


c.  Model  Alteration 

Model  alteration  capability  has  been  provided  that  allows 
the  user  to  adjust  arc  and  node  parameters  and  resolve  the  modified 
problem.  After  viewing  the  new  solution,  the  final  problem  attributes 
may  be  saved  to  a  data  file,  the  original  file  updated,  or  another  adjust 
ment  cycle  performed. 

d.  Reports 

The  program  displays  (on  user  designated  peripherals)  an 
echo  print  of  the  input  file  and  final  solution  reports.  Automatic 
generation  of  these  reports  is  the  default  case,  however,  either  report 
may  be  disabled.  Appendix  B  contains  examples  of  these  reports. 


V.  COMPUTATIONAL  RESULTS 


This  section  describes  the  capabilities  of  the  package  with  respect 
to  problem  size,  solution  speed,  and  versatility.  As  with  most  applied 
mathematical  programming  projects,  availability  of  suitable  example 
problems  hindered  the  testing  effort.  Although  standard  test  problems 
have  been  described  in  the  literature,  most  are  either  very  small  academic 
examples  or  randomly-generated  problems.  The  difficulty  with  using 
randomly-generated  problems  is  that  they  are  unstructured  and  therefore 
may  not  adequately  exploit  the  efficiency  of  a  programming  code  designed 
to  solve  "real-world"  problems.  It  has  been  suggested  [Ref.  3]  that  such 
random  problems  may  even  be  harder  to  solve  than  naturally  occurring 
formulations.  Nevertheless,- it  has  been  necessary  to  include  some 
randomly-generated  problems  to  provide  a  complete  assessment  of  package 
capabilities.  All  randomly  generated  problems  were  produced  by  NETGEN 
[Ref.  37]  running  on  an  IBM  370/3033  computer  and  subsequently  transferred 
to  the  microcomputer  using  package  file  transfer  features  (EDITOR  program 
and  commercial  data  transfer  program  [Ref.  35]).  Additionally,  a  few 
small  academic  examples  were  utilized.  Finally,  a  preprocessor  has  been 
written  to  construct  models  which  (although  fictitious)  more  closely  resemble 
real  world  situations  than  either  of  the  two  previously  mentioned 
problem  classes. 

A.  MAXIMUM  PROBLEM  SIZE 

The  maximum  representable  formulation  (in  terms  of  numbers  of  arcs 
and  nodes)  is  a  function  of  model  type,  the  relative  proportion  of  arcs 
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and  nodes,  and  the  microcomputer  employed.  All  statements  concerning 
problem  size  apply  to  the  APPLE  II  microcomputer  and  no  attempt  has  been 
made  to  translate  the  results  to  other  hardware  configurations. 

Although  one  may  specify  a  nominal  ratio  of  arcs  to  nodes  and  con¬ 
figure  the  programs  in  that  manner,  it  might  be  necessary  to  alter  this 
proportion  to  accommodate  a  particular  model.  A  limitation  of  this  package 
is  that  in  order  to  effect  such  an  alteration,  both  the  solution  program 
pertaining  to  the  class  of  models  in  question  and  the  MASTER  program  must 
be  recompiled.  The  coding  change  itself  is  trivial  (two  constants  in 
each  program  control  the  partition),  but  the  compilation  process  is  time 
consuming.  One  solution  to  this  dilemma  would  be  to  maintain  multiple 
copies  of  the  package,  each  with  different  problem  size  capabilities,  to 
represent  anticipated  modeling  requirements.  Table  1  gives  the  partition 
employed  during  the  development  of  the  package. 

The  APPLE  II  microcomputer  has  approximately  39,900  words  of  memory 
available  for  program  use  after  the  Pascal  operating  system  has  been 
loaded.  The  combined  code  of  the  main  segment  of  a  solution  program 
(always  in  memory)  and  that  of  its  largest  subordinate  segment  (typically 
the  solution  segment)  will  determine  the  memory  available  for  data 
storage.  In  the  absence  of  data  partitioning,  this  also  defines  the 
maximum  size  of  a  representable  model.  The  memory  budget  for  each 
solution  program  is  shown  in  Table  2. 

Data  storage  memory  requirements  can  be  expressed  as  a  function  of 
the  number  of  arcs  and  nodes  represented  and  will,  of  course,  differ  for 
each  solution  program.  These  relationships  are  as  follows: 
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TABLE  1 


NOMINAL  PACKAGE  DATA  PARTITION 


Nodes 

Arcs 

GNET 

400 

2000 

ELASTICNET 

300 

500 

NLPNET 

50 

175 

TABLE  2 

PACKAGE  MEMORY  USAGE 


Main 

Largest 

Auxiliary 

Available 

Segment 

Segment 

for  Data 

GNET 

5072 

7517 

27311 

ELASTICNET 

6565 

13898 

19442 

NLPNET 

6804 

19941 

13155 

Requirements  are  given  in  BYTES  (8-bits)  with  39,900 
total  BYTES  available  for  program  use. 
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1. 

GNET: 

6 

1  A  | 

1  ♦ 

20  | 

N  1 

£ 

27311, 

2. 

ELASTICNET: 

14 

1  A  | 

1  + 

26  | 

N  1 

< 

19442, 

3. 

NLPNET: 

56 

1  A  I 

1  + 

22  | 

N  1 

< 

13155, 

where  |  A  | 

signifies  the 

number 

of 

arcs 

and 

1 

N  |  the  number  of  nodes 

represented.  A  summary  of  this  information  for  selected  combinations  of 
arcs  and  nodes  is  presented  in  Table  3. 

B.  SOLUTION  TIMES 

Solution  times  vary  with  the  complexity  of  the  model  under  considera¬ 
tion.  The  simplest  formulation,  inelastic  linear  models  (GNET),  requires 
only  integer  arithmetic  for  logical  comparisions  and  therefore  produces 
the  fastest  solution  times  for  a  given  model  size.  Elastic  models, 
although  able  to  take  advantage  of  some  of  the  efficiencies  associated 
with  linear  models,  represent  a  versatile  but  complex  formulation  that 
necessitates  additional  work  to  cope  with  the  increased  information 
requirements.  By  specifying  "1-u"  arcs  in  a  model ,  numerous  subproblems 
(each  equivalent  to  a  single  ENET  run)  must  be  solved  and  coordinated. 

The  maximum  number  of  such  enumerations  (worst  case)  grows  exponentially 

L  ^  1 

with  the  number  of  "1-u"  arcs  (2  +  1,  where  k  is  the  number  of 

"1-u"  arcs),  hence  solution  time  increases  proportionally.  Nonlinear  models 
require  the  most  time  to  achieve  optimality.  These  models  employ  vast 
amounts  of  real  arithmetic  (including  transcendental  function  computations 
which  are  very  time  consuming  on  the  APPLE  II)  to  perform  the  necessary 
functional  evaluations,  direction  finding  and  linesearch  solutions. 

Nonlinear  solution  times  are  also  highly  sensitive  to  the  form  of  the 
objective  functions. 
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TABLE  3 


ALLOWABLE  ARC/NODE  PARTITIONS 


ARCS 


NODES 

GNET 

ELASTICNET 

NLPNET 

20 

3368 

1351 

227 

30 

3346 

1333 

223 

40 

3323 

1314 

219 

50 

3301 

1295 

215 

60 

3278 

1277 

211 

70 

3256 

1258 

207 

80 

3233 

1240 

203 

90 

3211 

1221 

199 

100 

3188 

1203 

185 

150 

3076 

1110 

175 

200 

2963 

1017 

250 

2851 

924 

300 

2738 

831 

350 

2626 

738 

400 

2513 

645 

450 

2401 

553 

500 

2288 

550 

2176 

600 

2063 

650 

1951 

700 

1838 

750 

1726 

800 

1613 

850 

1501 

900 

1388 

950 

1287 

1000 

1163 
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A  wide  variety  of  test  problems  have  been  examined  and  the  solution 


results  are  presented  in  Table  4.  The  solution  times  are  all  quite 
reasonable  (although  orders  of  magnitude  slower  than  those  obtained  with 
large  computers)  and  reflect  the  expected  relationships  between  model 
classes.  An  example  of  the  accuracy  for  a  representative  nonlinear 
formulation  is  given  in  Table  5. 

C.  MOOELING  FLEXIBILITY 

To  demonstrate  the  flexibility  of  the  package,  a  preprocessor  program 
has  been  written  to  model  a  three  echelon  product  distribution  network. 

In  such  a  model,  products  flow  from  production  centers  to  customers  via 
distribution  centers  with  storage  or  handling  costs  incurred  at  the 
intermediate  echelon.  The  object  of  this  formulation  is  to  meet  demands 
from  available  supplies  at  minimum  transportation  and  handling  cost. 
Transportation  costs  are  modeled  by  arcs  between  components  of  the 
various  model  echelons.  Handling  costs  can  be  depicted  by  adding  an 
additional  arc  at  each  distribution  center  with  appropriate  cost  function 
and  flow  bounds.  Figure  33  shows  the  basic  structure  of  such  a  model. 

The  preprocessor  requires  geographical  coordinates  and  product  flow 
requirements  for  each  location  to  be  modeled.  All  the  transportation  arcs 
are  assigned  linear  cost  functions  proportional  to  the  great  circle 
distance  between  the  various  locations.  Handling  arc  objective  functions 
are  assigned  by  the  user  to  reflect  the  desired  formulation.  The  user 
also  specifies  the  minimum  number  of  customers  to  be  linked  to  each 
distribution  center.  The  preprocessor  then  constructs  a  model  as  follows: 
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TABLE  5 

OPTIMAL  SOLUTION  COMPARISON 


Arc 


Optimal  Solution  Flow 


APPLE  II 


DEC  20/60 
[Double  Precision) 


1 

52.1063 

52.1063 

2 

47 .8936 

47.8937 

3 

18.5550 

18.5552 

4 

33.5514 

33.5511 

5 

26.4450 

26.4448 

6 

21.4486 

21.4489 

7 

5.000 

5.000 

8 

6.60537 

6.6050 

9 

6.9496 

5.9502 

10 

14.6698 

14.6613 

11 

10.0000 

10.0000 

12 

8.88158 

8.8898 

13 

0.0000 

0.0000 

14 

13.3946 

13.3950 

15 

13.0504 

13.0498 

16 

0.33020 

0.3387 

17 

0.0000 

0.0000 

18 

21.1184 

21.1102 

Optimal  Value  of 

Objective  Function 

1453.420 

1453.417 

Objective  function  components  are  power,  linear,  and  hyperbolic  SIN 
functions.  (Problem  is  NLPNET3). 


i 


113 


Production 

Centers 


Distribution 

Centers 


Customers 


Fig.  33.  Three  Echelon  Distribution  Model  With  Handling  Costs 


STEP  1,  Production  centers  are  linked  with  each  distribution  center. 
STEP  2.  Handling  arcs  are  added. 

STEP  3.  Each  production  center  is  connected  with  the  required  number 
of  customers  selected  in  order  of  proximity. 

STEP  4.  Any  customer  not  linked  to  at  least  one  distribution  center 
by  STEP  3  is  connected  to  the  closest  distribution  center. 
The  following  types  of  handling  cost  functions  can  be  accommodated: 

1.  Linear. 

2.  "1-u". 

3.  Fixed  cost. 

4.  Nonlinear. 

The  linear  formulation  results,  of  course,  in  a  completely  linear 
model.  Addition  of  "1-u"  arcs  transforms  the  model  into  one  where 
distribution  centers  are  either  open  at  full  capacity  (with  linear 


'jk*m »£/*•- ;< 
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handling  costs)  or  closed.  Fixed  cost  formulations  incur  an  additional 

charge  just  for  opening  a  distribution  center  (fixed  cost)  with  a  linear 

handling  cost  (variable  cost)  applied  to  each  unit  of  flow.  Such  fixed 

cost  models  are  generated  using  "1-u"  arcs  as  shown  in  Figure  34. 

Finally,  the  handling  costs  may  be  represented  by  various  nonlinear  cost 

functions  (potentially  different  for  each  handling  arc).  Two  nonlinear 

functions  that  immediately  come  to  mind  as  appropriate  in  a  handling  cost 

situation  are  quadratic  functions: 

Cost  *  f(x)  *  ax^  +  bx  +c 

lower  bound  <  x  <  upper  bound, 

and  hybrid  linear/quadratic  functions: 

cx  lowerbound  <  x  £  changeover  point. 

Cost  »  f ( x )  =  2 

ax*  +  bx  +  (c) (changeover  point) 

changeover  point  <  x  <  upperbound, 
a  ^  0  (for  convexity) . 
x  s  flow 

For  demonstration  purposes,  a  series  of  problems  were  generated  with 
the  following  structural  characteristics: 

1.  five  production  centers, 

2.  five  distribution  centers, 

3.  twenty-five  customers,  and 

4.  a  minimum  of  eight  customers  per  distribution  center. 

The  resulting  40-node,  73-arc  model  was  replicated  using  the  various 
handling  arc  objective  functions  described  above.  Solution  times  for 
these  test  problems  are  given  in  Table  6. 

Numerous  modifications  of  the  models  are  possible  by  exercising  the 
"elastic"  features  of  ENET  and  UNET.  In  this  manner,  realistic  enhancements 
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fixed  cost  ♦  variable  cost  x  ( upperbound -lowerboundj 
f(x)  *  ““  upperbound 


{0,  Upperbound} 


f(x)  *  -  Variable  cost 


0,  upperbound- lowerbound 


f  i  xed  h 
cost  I 


uppert 


Fig.  34.  Modeling  Fixed  Costs 

such  as  optional  surplus  or  shortage  assignment  can  easily  be  incorporated 
into  the  model. 


TABLE  6 


THREE  ECHELON  DISTRIBUTION  MODEL  SOLUTION  TIMES 
(40  Nodes,  73  Arcs  With  5  Handling  Arcs) 


Form  of  Handling  Arc 

Objective  Function 

Solution 

Time  (Seconds) 

Pivots  or 
(Iterations) 

Linear  (ENET) 

40 

92 

One  "1-u" 

120 

256 

Two  "1-u" 

175 

362 

Five  "1-u" 

584 

1179 

One  Fixed  Cost 

115 

260 

Two  Fixed  Cost 

180 

387 

Five  Fixed  Cost 

458 

960 

Quadratic 

65 

(4) 

Hybrid  Linear /Quadratic 

68 

(5) 

-  i 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  research  effort  described  by  this  paper  has  been  quite  successful. 
Microcomputer  adaptations  of  advanced  algorithms  for  minimum  cost  network 
flow  problems  have  been  shown  to  be  not  only  feasible  but  also  practical 
for  moderate-sized  formulations.  Additionally,  an  integrated  repertoire 
of  solution  methods  has  been  presented  that  illustrates  the  usefulness  of 
such  packages  and  serves  as  a  prototype  for  their  implementation  on 
larger  computer  systems.  This  software  fills  a  void  in  the  existing 
spectrum  of  computer  techniques  for  network  problems  by  providing  a 
single  package  to  deal  with  a  variety  of  diverse  modeling  situations. 

Economic  justification  for  the  use  of  such  microcomputer-based 
packages  requires  further  investigation  and  could  easily  be  the  exclusive 
subject  of  a  formal  study.  Certainly,  the  microcomputer  will  not  replace 
larger  computers,  but  instead  will  serve  as  a  useful  supplement.  The 
capabilities  of  this  software  package  infer  that  the  microcomputer's 
niche  lies  in  providing  quick  answers  to  relatively  small  problems.  In 
this  sense,  a  desktop  computational  device  might  be  more  convenient  to 
use  than  the  services  of  a  large  computer  center.  This  is  especially 
true  if  the  software  is  user-friendly  and  easy  to  use. 

There  are  certain  situations  where  the  microcomputer  is  clearly 
the  only  choice.  Consider  remote  sites  where  access  to  large  computers 
is  limited  or  nonexistent.  In  such  a  scenario,  the  portability  of  the 
smaller  computer  provides  the  analyst  with  a  powerful  mathematical 
programming  capability  that  would  otherwise  not  be  available.  Indeed, 


118 


microcomputers  have  already  been  used  in  primitive  locations  utilizing 
rudimentary  power  supplies,  so  such  a  proposal  is  not  idle  conjecture. 

The  finite  time  horizons  imposed  on  this  project  necessarily  limited 
the  scope  of  the  initial  study  and  there  are  many  enhancements  that  could 
be  pursued.  First,  the  algorithms  presented,  although  quite  efficient, 
could  be  improved  and  further  tuned  for  the  microcomputer  environment. 

Also,  the  use  of  partitioned  data  structures  was  not  investigated  as  a 
means  to  solve  larger  problems  with  available  resources.  Successful  use 
of  such  procedures  could  extend  the  usefulness  of  the  package.  Finally, 
the  most  obvious  extension  would  be  the  inclusion  of  additional  algorithms 
to  accommodate  other  classes  of  network  models.  For  example,  a  generalized 
network  code  recently  made  available  to  us  by  McBride  [Ref.  38]  would  be 
a  perfect  complement  to  the  algorithms  already  included. 

It  is  hoped  that  the  success  of  the  work  presented  here  will  further 
stimulate  the  development  of  additional  mathematical  programming  software 
for  microcomputer  use.  As  smaller  and  more  capable  computers  are  inevi¬ 
table,  the  operations  research  community  must  be  prepared  to  exploit 
their  considerable  potential. 
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APPENDIX  A 
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PACKAGE  EXTERNAL  DATA  STRUCTURE  (DATA  FILES) 


CONST 

MXNUMCOEF  =  3; 
MXNAMELENGTH  =  9; 


TYPE 

SMALLSTRING 

COEFARRAY 

RECTYPE 

NETWORK 

NODETYPE 

NETRECORD 


STRING[MXNAMELENGTH] ; 

ARRAY[1.. MXNUMCOEF]  OF  REAL: 

(HEADER  ARC  NODE)* 

( LINEAR ^NONLINEAR | ELASTIC, MIXED  INTEGER, GENERAL) ; 
( SUPPLY ,OEMAND , TRANSSHIPMENT) ; 

RECORD 

CASE  ENTITY:  RECTYPE  OF 
HEADER:  (PROBLEMNAME  :  SMALLSTRING; 

PROBLEMTYPE  :  NETWORK; 

NUMNODES  :  INTEGER; 

NUMARCS  :  INTEGER; 

DATECREATED  :  SMALLSTRING; 

DATELASTUPDATE  :  SMALLSTRING); 


ARC: 


(ARCNAME 

SOURCENODE 

DESTINATIONNODE 

LOWERBOUND 

UPPERBOUND 

INITIALFLOW 


SMALLSTRING; 

INTEGER; 

INTEGER; 

REAL; 

REAL; 

REAL; 


CASE  ARCTYPE:  NETWORK  OF 


LINEAR 

NONLINEAR 


(UNITCOST 

ZEROUARC 

(TYPEFTN 

NUMCOEF 


COEF:  COEFARRAY); 


REAL; 

BOOLEAN); 

INTEGER; 

1.. MXNUMCOEF); 


NOOE: 


END 


(NODENAME 

NODENUMBER 

NODEKIND 

NETFLOW 

LOWERRANGE 

UPPERRANGE 

LOWERPENALTY 

UPPERPENALTY 


SMALLSTRING; 

INTEGER; 

NODETYPE; 

REAL; 

REAL; 

REAL; 

REAL; 

REAL); 


NETFILE  :  FILE  OF  NETRECORD: 


1 


APPENDIX  B 
TYPICAL  REPORTS 

This  appendix  contains  three  examples  illustrating  typical  reports 
generated  by  the  various  solution  programs.  These  examples  are: 

Example  1:  Linearization  of  a  small  nonlinear  formulation  with  an 
infeasible  starting  solution.  This  serves  as  an  example  of  both  MASTER 
program  and  GNET  output. 

Example  2:  Typical  NLPNET  run  of  a  feasible  (starting  solution) 
nonlinear  formulation. 

Example  3:  Output  from  a  small  linear  formulation  solved  with  ELASTICNET. 
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DATE:  6  SEP  81 


APPL?-NET_:  SOLUTION  NODOSE 


Version  IIA 


Aug  8 1 


PROCESSING  DATA: NLPBETl .  NET  »  A  BON  LINEAR  NETWORK  PROGRAB. 
BOBBER  OP  NOOES  »  5  BOBBER  OP  ARCS  *  8 

PROBLEB  IS  INFEASIBLE 
LINEARIZATION  PERPORHED: 

7.  Assignment  or  zero  costs  for  all  arcs. 

THETA  (POINT)  *  0. 00000 


SO LOTION  NODDLE  CRAIN  OP  EVENTS: 
GNBT-->  NLPNET — >  SOLOTION. . . 


APPLEHET  -  GNET  BODOLE 
VERSION  IIIA  OP  28  A0G  81 


2000  ABC,  400  NODE  VERSION 

PILE:  DATA: NLPNET1.NET  CREATED:  1-SBP-81 
PROBLEB  BARE  IS  NLPNET- 1 

BOBBER  OP  NODES  «  5  BOBBER  OP  ARCS  *  8 


OPDATED:  1-SRP-81 


ARC  LIST. 

ARC 

NABE 

ONE 

TNO 

THREE 

POOR 

FIVE 

SIX 

SEVEN 

EIGHT 


PROB 

NODE 

1 

1 

2 

2 

2 

3 

3 

5 


TO 

NODE 

2 

3 

3 

4 

5 

4 

5 
4 


OBIT 

COST 

0.00 

8:88 

8:88 

8:88 

0.00 


LONER 

BOOND 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 


OPPER 

BOOND 

80.00 

60.00 

50.00 

tt5°o:0o°o 

SO'.OO 

60.00 


NODE  LIST  . 

•  « 

FOR  THOSE  NODES  EXPLICITLX  IN 

NODE 

NODE 

NET 

NODE 

NABE 

BOBBER 

PLOW 

TTPE 

S0P1 

1 

100.00 

SOP PLI 

ORAHED 

2 

0.00 

TRANS SHIP BENT 

ON ABED 

3 

0.00 

TRANSSHIPBENT 

DEB  1 

4 

-30.00 

DERAND 

DEB2 

5 

-70.00 

DEBAND 

SOPPLIES  > 

100 

DEHAND  »  100 

INITIAL 

0.00 


7  Pivots  perforsed 

PLOW 


TO 

NODE 


PROB 

NODE 

1 

2 

l 


OPTIBAL  COST 


PINAL  SOLOTION 


F(X  ) 


n 

30 

31 

0.00 


°0 

0 

8 
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